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Abstract 


In  the  first  year  of  a  planned  three  year  program  the  research  hats  aimed  at  defining 
data  structures  and  completing  tools  for  a  new  flexible  approach  to  scientific  programming 
and  problem  solving.  Both  problems  of  program  complexity  associated  with  changing 
models  and  physics  as  well  as  joined  and  disjoint  multiple  independent  patched  domain 
decompositions  for  treating  complex  geometries  can  be  systematically  organized  within  the 
context  of  the  directed  graph  programming  concept  being  explored.  Problems  and  parts 
of  problems  having  geometric  connectivity  or  its  analogs  such  as  precedence  relationships 
are  naturally  exhibited  and  easily  debugged  or  modified  in  the  graph.  The  solution  of 
problems  is  literally  to  traverse  the  graph.  For  the  planned  prototype  aerodynamic  simu¬ 
lation  facility  the  graphs  which  are  to  embody  grid  generation  and  Navier-Stokes  solution 
procedures  are  to  be  constructed  with  a  graphical  editor  hosted  in  a  high  performance 
graphics  workstation. 


The  efficient  global  data  structure  for  the  system  is  a  set  of  large  linear  arrays  in 
which  the  data  for  the  independent  quadrilateral  blocks  of  mesh  are  sequentially  stacked. 
The  directed  graph  is  to  control  procedures  that  fill  pointer  arrays  to  the  data  structure. 
We  implement  both  upwind  characteristics  based  implicit  computed  boundary  point  pro¬ 
cedures  for  physical  and  external  domain  boundaries  and  on  internal  patch  boundaries 
frozen  boundary  data  that  is  obtained  from  the  solution  on  adjoining  overlapping  meshes. 
The  approach  admits  segmented  grids  with  steps  and  holes  and  boundary  conditions  that 
are  nonhomogeneous  by  type. 


Additions  and  improvements  to  our  upwind  computational  tools  have  included  the 
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1.  Introduction 


As  major  strides  have  been  made  in  the  past  several  years  in  the  speed  and  memory 
of  computers,  similar  advances  have  been  achieved  in  details  of  computational  technique 
for  constructing  finite  difference  meshes  and  for  solving  the  compressible  Euler  or  Navier- 
Stokes  equations.  The  combined  progress  has  been  such  that  the  CFD  community  is  now 
challenged  by  the  possibility  of  solving  whole  interacting  flows  of  multicomponent  systems. 
Topical  examples  are  flows  about  aircraft  and  missiles,  and  in  turbomachines  and  rocket 
engines. 

Such  thinking  is  symptomatic  of  the  fact  that  the  research  focus  is  rapidly  shifting 
from  2-D  explorations  and  demonstrations  to  viable  3-D  computational  methodologies. 

The  shift  poses  two  research  challenges.  The  first  is  to  bring  a  new  generation  of 
algorithms  to  an  effective  fruition  in  terms  of  accuracy  and  computational  efficiency.  The 
second  is  to  organize  the  problems  of  constructing  meshes  and  solving  the  gasdynamic 
PDE’s  on  geometrically  complex  flow  domains  in  such  a  way  as  to  minimize  the  work  of 
problem  setup  and  execution  for  each  new  geometry  or  mesh  structure.  This  requires  the 
development  of  computational  systems  to  unify  the  many  components  of  the  simulation  - 
principally  grid  generation,  computation  and  graphics. 

The  present  program  contains  both  the  above  aspects  in  research  aimed  at  substan¬ 
tially  enhancing  accuracy  and  productivity  in  the  numerical  simulation  of  aerodynamic 
flows  in  complex  geometry.  The  program  builds  on  and  extends  recent  work  in  algebraic 
grid  generation1-3  for  composite  patched  mesh  systems  and  the  CSCM  upwind  implicit 
algorithm4-9  for  solving  the  compressible  Euler  and  Navier-Stokes  equations.  With  the 
conservative  upwind  method  we  have  successfully  demonstrated10  the  capability  to  ac¬ 
curately  capture  weak  reflecting  shocks  in  compressible  flow  using  solution  adaptive  grid 
refinement  strategies  with  aligned  overset  mesh  patches  like  those  proposed  by  Berger  and 
Oliger11.  On  conventional  grids,  the  implicit  upwind  method  exhibits  robust  stability 
properties  with  well  posed  implicit  characteristic  procedures4  at  all  external  flow  bound¬ 
aries.  However,  the  upwind  method  particularly  excells  in  accuracy  and  stability  in  the 
vicinity  of  interior  overset  patch  boundaries  with  only  a  minimal  overlap  as  required  to 
exchange  data  through  interpolation6. 

The  combined  strategy  of  employing  upwind  methods  with  both  composite  and  over¬ 
set  patched  meshes  is  sufficient  to  treat  problems  of  any  geometric  and  flow  structure 
complexity  with  desired  accuracy  and  computational  efficiency. 

Mesh  generators  and  flow  solvers  can  be  programmed11,12  in  FORTRAN  with  a  flexi¬ 
ble  global  data  structure  comprising  a  linear  array  with  pointers  and  parameter  controlled 
boundary  interface  routines  to  handle  general  multiple  mesh  decompositions  of  the  domain. 
However,  both  the  programming  of  the  codes  and  the  setup  to  solve  particular  complex 
problems  is  tedious  and  invites  mistakes  that  may  be  hard  to  track  down.  The  problems 
are  inherent  in  the  unnatural  structure  of  linear  programming  languages  like  FORTRAN 
that  do  not  exhibit  the  connectivity  of  objects. 
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Recent  research  by  Oliger  and  his  students  has  supported  the  possibility  to  write 
programs  and  execute  them  in  a  more  natural  way  that  directly  represents  the  connectivity 
of  objects.  The  approach  is  based  on  the  directed  graph  which  is  to  be  implemented  and 
exhibited  in  a  workstation  based  graphical  interface.  The  graph  is  a  natural  representation 
of  the  macro  data  flow  of  the  program.  An  important  feature  of  the  system  is  that  it  is 
recursive,  i.e.  nodes  of  the  graph  may  represent  pieces  of  the  graph.  To  facilitate  using 
programs  with  disparate  data  structures,  filters  are  inserted  between  nodes  to  reformat  the 
output  (range)  data  of  one  process  to  be  suitable  for  the  input  (domain)  of  the  next.  A 
language  such  as  LIL,  proposed  by  Goguen13,  is  being  developed  for  writing  these  filters. 
To  minimize  programming  while  permitting  great  flexibility,  programming  will  be  carried 
out  and  programs  executed  in  the  context  of  the  UNIX  environment.  The  Make  facility  of 
UNIX  naturally  lends  itself  to  achieve  the  composition  of  these  programs. 

Based  in  the  above  technologies,  a  major  purpose  of  the  evolving  program  will  be 
to  construct  and  test  a  prototype  programming  system  for  aerodynamic  simulation.  The 
simulation  will  generate  multiple  patch  meshes  for  representative  problems  and  direct  the 
flow  solver  over  the  meshes. 

Other  complementary  facets  of  the  research  are  to  seek  further  improvements  in  ac¬ 
curacy  and  efficiency  of  the  upwind  method  and  algebraic  grid  generation  procedures  that 
we  have  been  studying. 

Topics  that  we  are  pursuing  in  the  upwind  methods  area  that  derive  from  our  recent 
research  are,  first,  strongly  conservative  schemes  with  second  order  spatial  accuracy  in 
the  complete  2-D  and  3-D  Taylor  series’  expansion.  Here,  second,  these  schemes  are 
expressed  in  available  families  of  upwind  biased  finite  volume  computational  cells  rather 
than  the  classical  data  node  centered  cells.  On  such  cells  an  interesting  hierarchy  of 
splittings  and  biased  upwind  methods  are  possible  that  are  related  to  and  extend  the 
locally  one  dimensional  higher  order  TVD  schemes  of  Chakravarthy  and  Osher14  and 
Harten15.  Related  flux  limiting  schemes  including  some  based  on  the  work  of  Yang,  et  al.8 
are  to  be  explored. 

Another  research  topic  we  are  embracing  and  one  that  has  recently  shown  great 
promise,  particularly  for  steady  flow  implementation  with  upwind  schemes9,  is  flow  struc¬ 
ture  adaptive  mesh  redistribution.  This  technique  requires  no  additional  patches  of  mesh 
(with  additional  computational  costs)  to  achieve  the  same  order  of  magnitude  improvement 
in  resolution10  as  is  achieved  with  mesh  embedding  and  overset  refinement. 

Topics  in  the  algebraic  grid  generation  area  emphasize  extension  to  3-D.  One  impor¬ 
tant  facet  is  the  construction  of  smooth  body  conforming  surface  patches  of  curvilinear 
coordinate  mesh.  Our  present  mesh  generation  philosophy,  which  is  suitable  for  both  2-D 
and  3-D,  emphasizes  techniques  based  on  low  order  polynomial  blending  of  parametric 
cubic  and  tension  spline  functions.  These  are  used  for  both  bounding  and  coordinate 
curves.  Blended  stretching  functions  are  used  to  obtain  desired  mesh  point  distributions. 
The  approach2  represents  a  simplification  of,  and  borrows  many  tools  from,  the  transfi- 
nite  interpolation  method  of  Vinokur  and  Lombard1.  The  latter  relates  to  the  work  of 
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Eriksson16  and  to  Eiseman  and  Smith17.  For  boundary  surface  meshes,  our  approach 
which  suppresses  the  explicit  dependence  on  boundary  gradients  of  Hermite  interpolation 
follows  the  work  of  Gordon18. 

Another  development  of  the  boundary  blending  philosophy  which  we  are  persu¬ 
ing  is  patched  mesh  generation  using  isoparametric  macro  finite  elements  (Baker  and 
Manhardt19).  The  basis  functions  for  the  elements  can  then  be  used  for  coordinate  trans¬ 
formation  (Mitchell  and  Wait20).  A  feature  of  the  macro  element  approach  which  has 
recently  been  advanced  in  2-D  by  Oliger  and  Suhr3  is  the  use  of  piecewise  quadratics  for 
boundary  curves  in  a  way  that  achieves  C1  continuity. 

An  important  topic  of  the  new  programming  system  that  we  are  presently  studying  is 
how  to  represent  the  algebraic  grid  generation  problem  for  composite  meshes  in  a  directed 
graph  in  such  a  way  that  the  connectivity  (continuity  properties)  are  elicited  and  the 
procedure  is  made  very  efficient.  A  similar  topic  exists  for  gasdynamics  in  how  to  efficiently 
direct  the  operationally  explicit  symmetric  Gauss-Seidel  implicit  solution  algorithm  on  the 
composite  mesh. 

Finally,  as  suggested  above,  computer  graphics  naturally  fits  within  the  context  of 
a  complete  simulation  system.  From  an  earlier  effort  we  have  a  2-D  graphics  package 
programmed  in  a  flexible  multiple  mesh  data  structure  to  mirror  that  of  the  solution 
procedure12. 

In  the  following  sections  we  provide  additional  details  about  the  implementation  of  the 
directed  graph  based  programming  environment  and  advances  in  computational  methods 
that  we  have  made. 
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2.  Directed  Graph  Programming  Environment 

Much,  if  not  most,  effort  in  science  is  devoted  to  communication  between  various  media 
and  processes.  These  links  are  man-machine,  machine-machine,  and  program-program, 
etc.  The  purpose  of  this  project  is  to  develop  an  extensible  framework  in  which  this 
communication  can  efficiently  take  place  and  these  processes  controlled. 

Man-machine  interaction  is  especially  costly  and  tedious  when  geometrically  compli¬ 
cated  problem  domains  need  to  be  specified.  The  use  of  new  architectures  and  algorithms 
or  the  introduction  of  new  physics  often  pose  severe  time  and  cost  penalties  because  of 
extensive  software  modifications.  These  two  problems  are  being  focused  upon  initially  and 
the  system  will  then  be  expanded  in  its  role  and  functionality. 

The  system  is  envisioned  to  encompass: 

1.  Graphical  input  via  a  graph  structured  terminal  interface 

2.  Graph  structure  top  down/bottom  up  development  system  which  facilitates  combining 

existing  code  in  different  languages 

3.  An  object  oriented  interprocess  communication  system  with  several  levels  of  formality 

4.  Formal  manipulation  ability  for  derivations 

Program  Representation,  Macrodataflow  and 
Directed  Graphics 

The  major  obstacle  to  building  large  programs  from  program  fragments,  and  perhaps 
utilizing  pre-existing  programs,  is  the  management  of  data  and  data  structures.  The 
program  itself  will  have  defined  input  and  output  and  the  program  can  be  visualized  as 
the  process  of  performing  a  sequence  of  functions  on  the  input  in  order  to  obtain  the  desired 
output.  A  program  can  be  organized  from  the  top  down  by  specifying  these  functions  and 
their  relationship  to  one  another.  If  the  program  is  broken  into  rather  large  pieces  called 
nodes  involving  several  operations  (functions)  on  several  operands  (data),  then  we  can 
describe  this  process  in  terms  of  macro  data  flow.  For  each  node  we  describe  the  data 
required  for  each  functional  operation,  the  functions  to  be  performed,  and  the  output  to 
be  produced.  We  then  complete  the  specification  of  the  program  by  linking  these  nodes 
together  in  the  required  order.  This  process  can  be  described  nicely  using  the  idea  of 
directed  graphs.  A  directed  graph  consists  of  a  collection  of  nodes  N  =  {n  :}  and  links  or 
edges  E  =  {e<y}.  The  edges  =  (t, j)  are  ordered  pairs  with  tail  »  and  head  j  leading 
from  node  n*  to  node  ny  and  are  unidirectional  paths.  The  matrix  E  =  {etJ}  made  up 
of  the  array  of  edges  is  not  generally  symmetric.  The  nodes  can  be  used  to  define 
the  processes  and  the  edges  used  to  describe  the  flow  of  data  through  the  system.  In  a 
labeled  directed  graph  labels  /,•  are  attached  to  nodes  to  describe  their  function.  On  a  very 
primitive  level  we  can  describe  the  computation  of  d  =  (o  4-  b)  •  c  as 
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The  first  major  step  in  the  construction  of  a  modern  compiler  is  the  construction  of  a 
directed  graph  which  represents  the  computation.  We  plan  to  graphically  describe  the 
functions  and  their  relationships  to  each  other  on  a  traditional  subroutine  or  subprogram 
modularity.  The  graph  then  needs  no  aritificial  linear  ordering.  Its  struction  will  be 
coupled  with  an  abstract  data  description  language  which  will  provide  a  facility  for  the 
generation  of  filters.  These  filters  will  take  output  from  existing  programs  and  modify  it 
for  input  to  the  succeeding  node  as  required.1 

Visual  Graph  Generating  User  Interface 

We  have  well  in  development  a  visual  communication  system  based  upon  windows  and 
menus  such  as  those  used  on  the  SUN  workstations.  We  intend  to  extend  this  capability 
by  implementing  a  heirarchy  of  menus  in  a  tree  structure.  In  the  program  development 
environment  the  relationships  of  processes  are  explicitly  defined  by  connecting  windows  to 
one  another  via  a  directed  graph  which  then  serves  as  the  (executable)  program  control. 
No  linear  program  is  required  and  this  organization  can  be  done  differently  for  different 
hardware  systems  to  take  advantage  of  parallel  or  vector  processing,  etc.  The  basic  opera¬ 
tions  are  the  construction  of  a  graph  via  windowing  software,  the  “execution”  or  traversal 
of  such  a  graph,  and  the  graphical  display  of  the  graph  for  the  user  during  definition 
and  modification.  This  is  an  object  oriented  system  with  automatic  compatibility  of  do¬ 
mains  and  ranges,  etc.  Each  “window”  may  contain  subwindows  or  programs  in  different 
languages,  initially  at  least  FORTRAN,  PASCAL  and  C. 

Interprocess  Communication 

Communication  between  logical  sections  of  programs  (windows)  will  be  object  oriented 
and  file  based.  However,  several  levels  of  formality  will  be  available.  An  informal  level  like 
that  afforded  by  UNIX  pipes  will  be  followed  by  more  formal  files  with  internally  defined 
structures  and  features  which  allow  traces,  etc.  Finally,  formal  communication  with  a 
permanent  database  will  be  implemented  and  be  message  based.  Our  implementation  will 
be  a  macro  dataflow  system  without  interprocess  side  effects. 

Our  more  formal  files  with  internally  defined  structures  ;vill  be  used  to  communicate 
and  translate  data  between  formats  and/or  data  structures  to  allow  the  interaction  of 
programs.  This  allows  us  to  provide  an  upward  extensible  but  downward  compatible 
environment  when  we  enable  processes  to  read  files  extracting  only  those  segments  that 
they  “understand”. 
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Implementation 

We  plan  to  implement  this  system  in  the  spirit  in  which  it  is  planned  to  operate  ~ 
it  will  be  built  up  using  as  much  available  software  as  possible.  Pieces  will  be  redesigned 
only  when  it  becomes  obvious  that  it  is  desirable  to  do  so. 

Our  system  will  be  UNIX  based  and  interact  with  programs  at  the  operating  system 
level,  i.e.,  we  will  execute  our  graphics,  not  merely  concatenate  them  as  files,  etc.  We  plan 
to  use  windowing  features  and  existing  workstations  and  software  to  build  from. 

In  particular,  we  plan  to  incorporate  the  MACSYMA  Symbolic  manipulation  system 
or  equivalent.  We  are  currently  searching  existing  sofware  for  graphics  and  designs  and 
evaluating  which  would  incorporate  best. 

References 

1.  A.  Aho,  J.  Hopcroft,  J.  Ullman:  Data  Structures  and  Algorithms ,  Addison-Wesley, 
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3.  Accurate  Numerical  Simulation  of  Supersonic  Jet  Exhaust  Flow  with 

CSCM  on  Adaptive  Overlapping  Grids 

Design  of  propulsive  systems,  drag  reduction  and  stability  in  re-entry  vehicles,  pro¬ 
jectiles,  missiles  and  other  supersonic  aerodynamic  bodies  require  design  and  analysis 
computational  tools  that  can  accurately  predict  the  overall  flow  and  in  particular  the  base 
flow.  The  numerous  numerical  schemes  available  today  have  their  own  virtues  and  lim¬ 
itations.  The  present  paper  outlines  and  applies  the  coupled  overlapping  grids  CSCM 
scheme1,2,3,4  and  self-adaptive  grid  scheme5,6  to  the  jet  exhaust  base  flow  problem.  The 
work  culminates  a  series4,7  of  numerical  experiments  in  the  effects  of  mesh  topology  and 
resolution  on  base  flow  prediction.  All  the  computational  experiments  have  been  made  with 
the  same  adaptation7  of  the  Baldwin-Lomax8  algebraic  eddy  viscosity  model  of  turbulent 
interaction.  Here  computed  adaptive  grid  solutions  are  presented  and  comparisons  with 
results  of  experiment  and  other  numerical  schemes  are  made.  Great  improvement  to  the 
flow  solutions  are  achieved  through  flow  structure  adaptive  regridding  and  the  computed 
base  pressure  distribution  compares  very  well  with  experimental  predictions  of  Petrie  and 
Walker9. 

In  recent  years  evolving  computational  techniques  have  been  used  to  predict  super¬ 
sonic  viscous  base  flows4,5,7-23  in  a  variety  of  situations.  Complexity  of  the  supersonic 
base  flow  features  and  the  topological  difficulty  of  the  flow  region  have  led  to  discretization 
of  the  flow  field  in  different  Ways4,7,9,11,15,17'19,23.  Most  numerical  schemes  through  their 
programming  have  been  limited  to  obtaining  solution  in  a  single  rectangular  computa¬ 
tional  domain.  Hence  the  base  flow  region  is  often  discretized  by  a  wraparound  grid  or 
with  rectangular  composite  block  grids.  Recent  advances  in  data  structure  and  boundary 
treatment  have  opened  the  way  to  novel  grid  techniques1,2,3,4,  such  as  multiple  patched 
and  overlapping  grids,  and  have  provided  alternatives  to  solving  the  flow  field  in  simple 
rectangular  computational  domains.  Finally,  the  use  of  auto-adaptive  grid  procedures  re¬ 
sults  in  appropriate  grid  point  redistribution  and  improves  solution  accuracy  greatly  with 
minimal  computer  and  human  effort.  The  resulting  improvement  in  local  flow  structure 
capture  in  turn  improves  the  overall  solution  accuracy. 

The  Tactical  Missile  Turbulent  Base  Flow  Problem 

A  cold  flow  model  tactical  missile  turbulent  base  flow  was  experimentally  studied 
at  AEDC  and  the  experimental  results  and  other  numerical  predictions  were  reported  in 
Ref.  9.  The  axisymmetric  flow  experiment  featured  a  tangent  ogive-cone-cylinder  body 
terminating  in  a  square  base.  The  propulsive  jet  flow  exhausted  into  the  base  region 
through  a  centered  transonic  nozzle  with  a  flush  exit.  The  nozzle  occupied  0.2  diameter 
at  the  exit  plane  and  the  model  was  9  diameters  long.  The  free  stream  Mach  number 
was  maintained  at  1.4  and  the  Reynolds  number  of  the  flow  based  on  base  diameter  was 
0.92  x  10°.  The  nozzle  jet  flow  at  the  exit  was  at  a  nominal  Mach  number  of  2.7  and 
the  nozzle  exit  pressure  ratio  was  2.15  relative  to  the  free  stream  pressure.  Other  relevant 
wind  tunnel  conditions  are  given  in  Ref.  7. 
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The  schematic  of  the  above  flow  is  shown  in  Figure  1.  The  supersonic  external  flow  and 
the  nozzle  exhaust  flow  both  undergo  strong  expansions  around  the  corner  regions  followed 
by  recompression.  The  flow  separation  in  the  base  region  requires  the  flow  to  turn  with 
a  recompression  shock  near  the  outer  flow  region  and  a  barrel  shock  near  the  axis  being 
formed.  In  addition,  strong  shear  layers  and  associated  counter  rotating  vortices  form  in 
the  base  region.  Mach  disk  and  reflected  shocks  also  occur  within  the  flow  domain.  The 
occurrance  and  interaction  of  the  above  mentioned  flow'  features  make  the  jet  exhaust  base 
flow  problem  very  interesting  but  difficult  to  compute.  As  is  well  known,  accurate  base 
pressure  prediction  has  been  very  difficult  to  achieve.  The  prediction  depends  very  much 
on  the  appropriateness  of  the  grids  and  the  accuracy  of  the  boundary  and  interior  difference 
schemes.  The  present  study  indicates  that  an  algebraic  model  of  turbulence7  with  roughly 
appropriate  length  scale  and  velocity  measures  can  perform,  perhaps  surprisingly,  w-ell. 

Computational  Approach 

Numerical  solutions  to  the  base  flow  problem  have  been  obtained  using  various  schemes 
and  Ref.  16  reviews  most  results.  Solutions  to  the  base  flow  problem  using  the  Conservative 
Supra  Characteristic  Method24  {CSCM)  with  a  wraparound  and  a  simple  Cartesian  step 
grid  are  presented  in  Ref.  7.  As  pointed  out  there,  it  is  not  possible  to  achieve  adequately 
balanced  resolution  with  simple  mesh  topologies  on  geometrically  complex  flow  domains. 

Overlapping  grids  are  natural  and  appropriate  for  complex  flows.  The  choice  of  locally 
topologically  appropriate,  multiple,  independent  overlapping  grids  alleviates  the  geometric 
difficulty3’4.  With  well  balanced  resolution,  the  resulting  global  solution  is  more  accurate. 
The  use  of  multiple  grids  simplifies  the  discretization  of  the  complex  physical  domain  and 
allows  the  easy  use  of  a  solution  adaptive  scheme.  Corner  regions  are  modelled  accurately 
with  appropriate  local  wraparound  grids. 

The  details  of  the  CSCM  scheme  for  multiple  overlapping  grids  and  the  interior  grid 
boundary  data  procedure  are  given  in  Ref.  1  and  2.  The  overlapping  grid  CSCM  data 
structure  takes  proper  care  of  information  transfer  between  grids.  The  stable  character¬ 
istics  based  boundary  data  transfer  procedure  and  the  interior  upwind  difference  scheme 
together  predict  smooth  solution  everywhere.  Shocks  and  other  discontinuities  pass  from 
one  grid  to  the  other  without  distortion.  Solutions  to  complex  shock  reflection  problems 
with  overlapping  grids  were  reported  in  Ref.  2  and  4. 

Once  a  converged  solution  is  obtained  on  the  initial  grids,  the  self  adaptive  grid  scheme 
of  Nakahashi  and  Deiwert5  is  applied  and  the  grid  points  are  moved  in  all  the  grids  ac¬ 
cording  to  the  specified  constraints  and  flow  variations.  The  adaptive  gridding  is  done 
in  all  coordinate  directions.  Then  starting  with  the  interpolated  solution  on  the  adapted 
grid,  converged  solutions  are  obtained  using  the  CSCM  procedure  once  again.  Further 
grid  adaptations  and  computation  cycles  are  done  if  necessary  but  not  employed  for  the 
results  presented  here. 
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Results  and  Discussion 


Seven  overlapping  grids  are  used  in  the  present  problem.  The  first  grid  covers  the 
transonic  nozzle  and  the  second  grid  extends  from  the  first  grid  down  towards  the  outflow 
boundary  and  covers  the  exhaust  jet.  The  third  grid  is  used  for  the  outer  flow  up  to  the  base 
and  the  fourth  grid,  similar  to  the  second  grid,  extends  from  the  third  grid  downstream 
to  the  outflow.  The  fifth  grid  covers  the  base  flow  region.  The  sixth  and  seventh  grids 
are  local  wraparound  corner  grids  at  the  lower  and  upper  base  corner  regions.  Figure  2 
shows  the  boundaries  of  each  grid  and  Figure  3  shows  the  overlapping  grids.  All  the  grids 
have  at  least  two  cel!  overlap  with  their  neighbors.  The  grid  points  in  grids  1,  3,  and  5  are 
clustered  to  the  wall  to  resolve  the  viscous  layer.  A  blow  up  of  the  grids  near  the  nozzle 
exit  plane,  given  in  Figure  4,  shows  the  grid  structure  in  the  vicinity  of  the  corner.  Each 
of  the  seven  grids  span  a  much  simpler  sub-region  of  the  global  physical  domain.  All  grids 
were  generated  using  an  interactive  algebraic  grid  generation  scheme25.  We  note  that  with 
the  technique,  since  each  sub-grid  spans  only  a  simpler  region,  it  would  be  very  easy  to 
implement  nozzle  or  external  design  (shape)  changes  without  changing  the  grid  globally. 

The  interface  boundary  conditions  for  the  upwind  schemes  have  been  solved  satis¬ 
factorily  and  are  addressed  in  Ref.  1.  In  this  paper  we  apply  the  characteristic  interpo¬ 
lated  boundary  treatment1  at  all  overlapping  interior  boundaries  and  appropriate  physical 
boundary  conditions24,26  are  applied  at  the  boundaries  of  the  computational  domain.  At 
all  the  wall  boundaries  adiabatic,  no  slip  boundary  conditions  are  imposed. 

The  flow  is  nearly  stagnant  at  the  inflow  to  the  nozzle  section  with  a  stagnation 
pressure  of  307.7  psi  and  stagnation  temperature  of  207.0°  F.  Experimentally  measured 
flow  conditions18  are  imposed  along  the  inflow  boundary  of  grid  3,  at  a  station  16  nozzle 
exit  radii  upstream  of  the  base.  The  outer  flow  Mach  number  is  1.343.  With  these  specified 
conditions  a  converged  steady  state  solution  is  obtained  using  the  CSCM  method.  In  the 
nozzle  and  in  the  outer  supersonic  regions,  thin  layer  Navier-Stokes  equations  are  solved. 
In  the  base  flow  region,  (grid  5)  viscous  stresses  along  both  computational  coordinate 
directions  are  included.  An  algebraic  turbulent  eddy-viscosity  model8  reported  in  Ref.  7 
is  adapted  here  to  model  the  turbulent  interaction.  In  grids  1,  3,  6  and  7  wall  turbulence 
model  is  used  and  in  grids  2,  4  and  5  a  wake  model  is  used.  A  solution  of  spatially 
second  order  accuracy  is  obtained  in  about  600  iterations.  Mach  contours  of  the  computed 
solution  without  mesh  adaptation  are  shown  in  Figure  5.  The  contour  plots  show  the 
strong  expansion  regions,  the  recompression  shock,  the  barrel  shock,  the  mach  disk,  the 
reflected  shock  and  the  shear  layers. 

Solution  Improvement  with  Adaptive  Grids 

Further  improvements  to  the  above  solution  without  the  addition  of  more  grid  points 
is  easily  achieved  with  the  self-adaptive  grid  technique  of  Nakahashi  and  Deiwert5.  The 
adaptive  grid  scheme  is  automatic,  simple  to  use  and  very  flexible.  Grids  can  be  adapted 
in  all  coordinate  directions  and  grid  adaptation  can  be  performed  globally  or  locally. 

The  self-adaptive  grid  technique  has  been  used  in  the  past  with  simple  and  block 
grids.  The  redistribution  of  grid  points  in  the  overlapping  grids  is  more  involved,  but  not 
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too  difficult  depending  on  the  choice  of  the  overlapping  grids.  With  adaptive  regridding, 
some  of  the  interior  patch  boundaries  (that  do  not  conform  to  physical  boundaries)  should 
deform  according  to  the  flow  features,  but  physical  boundaries  such  as  wall,  inflow,  axis 
and  outflow  boundaries  should  remain  fixed.  With  present  programming,  if  each  grid  is 
adapted  independently,  then  all  the  interior  boundaries  will  remain  fixed.  To  allow  the 
interior  boundaries  to  move,  a  block  of  grid  points  in  the  region  of  interest  are  assembled 
from  all  the  grids  and  these  grid  points  are  moved  according  to  the  constraints  specified. 
After  regridding,  a  reverse  transformation  which  preserves  the  original  overlap  is  performed 
and  individual,  overlapped,  adapted  grids  are  formed.  This  proceedure  does  not  alter  the 
basic  auto-adaptive  grid  code  or  the  flow  solver. 

Based  on  the  converged  solution  shown  in  Figure  5,  the  grids  after  many  global  and 
local  adaptations  from  the  grids  of  Figure  3  appear  as  shown  in  Figure  6.  The  adapted 
grids  show  all  the  prominent  features  of  the  flow  field  (Figure  5)  based  on  which  the 
adaptation  was  performed.  A  combination  of  Mach  number  and  density  variations  were 
used  in  developing  the  adapted  grids.  Since  the  adaptation  is  a  result  of  flow  variation 
and  user  specified  constraints,  care  must  be  taken  during  adaptation.  It  is  possible  to 
over-adapt.  A  rule  of  thumb  is  to  first  adapt  globally  and  then  adapt  locally  to  flow 
structures. 

Steady,  converged  solution  on  the  adapted  grid  is  obtained  in  another  200  steps.  Inter¬ 
polated  solution  from  the  non-adapted  grid  served  as  the  starting  solution.  The  converged 
adaptive  grid  solutions  are  presented  in  terms  of  pressure,  density,  Mach  and  tempera¬ 
ture  contours  in  Figures  7  through  10.  Comparing  Figures  5  and  9,  one  can  see  that  the 
adapted  grids  indeed  serve  to  capture  the  flow  better  and  the  flow  structure  details  are 
more  prominent  than  on  the  non-adapted  grid  solution.  Specially  the  shear  layer  features 
and  the  Mach  disc  and  the  reflected  shock  near  the  axis  are  captured  more  sharply.  The 
merging  of  the  shear  layers  is  best  seen  from  the  temperature  contour  plots.  Complex  flow 
interactions  between  the  shear  layers,  and  the  recompression  and  barral  shocks  determine 
the  size  and  shapes  of  the  separated  base  flow  near  the  annular  base  region.  The  predicted 
base  pressure  profiles  with  and  without  grid  adaptation  are  compared  with  the  experi¬ 
mentally  measured  pressures  in  Figure  11.  Base  pressure  predictions  from  other  numerical 
results  taken  from  Refs.  4  and  9  are  reproduced  here  in  Figure  12.  The  adapted  overlap¬ 
ping  grids  solution  predicts  base  pressure  significantly  better  compared  with  non-adapted 
solution  and  other  numerical  predictions.  As  is  well  known  ar.d  reinforced  by  the  present 
results,  the  base  pressure  is  very  sensitive  to  the  flow  details  in  the  base  regions.  It  is  our 
belief  that  the  adapted  grids,  particularly  through  accurate  capture  of  the  shear  layers, 
allows  the  flow  field  near  the  base  region  to  evolve  accurately  and  this  coupled  with  proper 
application  of  the  adiabatic  boundary  conditions  at  the  base  wail  recovers  the  correct  base 
pressure.  In  Figure  13,  the  velocity  vectors  are  shown.  Various  recirculating  vortices  are 
clearly  visible.  In  Figure  14,  stream  lines  as  particle  traces  are  shown.  The  particles  are 
released  all  along  the  exit  plane  and  the  resulting  particle  traces  give  further  insight  into 
the  physical  aspect  of  the  flow.  The  particle  released  from  the  top  half  of  the  annular 
regions  is  lifted  up,  and  turns  downward  and  reaches  the  nozzle  exit  corner.  From  there 
they  are  convected  downstream  through  various  vortices.  The  empty  regions,  representing 


the  stagnant  regions  show  the  location  of  slower  moving  vortices  region.  The  Mach  number 
profile  along  the  axis  from  the  non-adapted  and  adapted  solutions  are  shown  in  Figure  15. 
The  Mach  number  is  a  little  higher  than  1.0  behind  the  shock  reflection  point  (where  it 
should  be  1.0)  Through  further  local  refinemnet  and  adaptation,  the  Mach  disc  can  be 
captured  better. 

Concluding  Remarks 

The  self-adaptive  grid  scheme  of  Nakahashi  and  Deiwert  coupled  to  the  overlapping 
grids  CSCM  scheme  is  demonstrated  to  be  very  effective  in  obtaining  accurate  solutions 
to  the  complex  jet  exhaust  base  flow  problem.  The  adapted  overlapping  grid  solution 
captures  all  the  important  flow  features  accurately  and  the  flow  details  agree  very  well 
with  experimental  observation.  Further,  excellent  agreement  between  the  predicted  base 
pressure  and  the  experimental  results  increases  the  confidence  in  present  flow  solver  and 
grid  strategy.  The  flexibility  of  the  combined  adaptive  and  overlapping  grids  makes  the 
CSCM  flow  solver  a  very  effective  design  and  analysis  tool. 
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Figure  4.  Grids  in  the  vicinity  of  3  corner. 


Figure  6.  Adapted  overlapping  grids.  Figure  8.  Pressure  contours  from  the  adapted  grids 
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gure  14.  Stream  lines  in  the  base  region. 


4.  Accurate,  Staggard  Grid  Upwind  Schemes  on  Biased  Upwind  Cells 

If  one  sees  the  seventies  and  first  half  of  the  eighties  as  a  period  of  learning  how  to 
construct  stable,  accurate  and  computationally  fast  methods  to  solve  compressible  flow 
equations,  particularly  the  Navier-Stokes  equations,  then  it  is  quite  apparent  the  focus 
is  shifting  and  for  a  time  the  emphasis  will  be  on  synthesis,  application  and  enhanced 
productivity.  The  present  paper  introduces  improvements  in  the  multi-dimensional  for¬ 
mulation,  numerical  implementation  and  mode  of  application  of  the  implicit  Conservative 
Supra-Characteristics  Method  (CSCM).  The  method  for  both  the  compressible  Euler  and 
Navier-Stokes  equations  was  first  presented  for  inviscid  quasi  one  dimensional  flow  in  the 
8th  ICNMF  by  Lombard,  Oliger  and  Yang1.  Extension  of  CSCM  to  a  2-D  and  source 
term  free  axisymmetric  formulation  was  given  by  Lombard,  Bardina,  Venkatapathy  and 
Oliger2.  Reference  2  introduced  solution  of  the  difference  equations  by  a  two  data  level 
block  tridiagonal  procedure  based  on  a  diagonally  dominant  approximate  factorization 
DDADI.  Subsequent  extensions  of  that  approximate  factorization  were  shown  by  Lombard, 
Venkatapathy  and  Bardina3  to  form  the  basis  for  a  family  of  single  data  level  operationally 
explicit  implicit  relaxation  algorithms.  One  of  those  algorithms,  a  symmetric  Gauss-Seidel 
space  marching  procedure  implemented  as  a  method  of  lines  in  two  space  dimensions  (or 
planes4  in  three  dimensions)  was  found  to  be  an  order  of  magnitude  more  rapidly  conver¬ 
gent  to  steady  state  than  the  same  equations  solved  in  the  two  level  linearized  pseudo  time 
relaxation  procedure.  In  the  generalized  Roe  form5 

AAq=AF  (1) 

which  relation  is  preserved  by  the^CSCM  constructed  matrix,  the  adjacent  grid  point 
interval  averaged  Jacobian  matrix  A  onto  the  conservative  variable  difference  is  used  to 
represent  the  flux  difference.  Thus  updating  A  and  q  in  the  locally  iterated  pseudo  time 
relaxation  procedure  at  each  marching  step  is  equivalent  to  updating  the  flux. 

Heretofore  we  wrote3  the  method  (here  sketched  as  a  first  order  single  level  method 
in  one  dimension)  as 

D  6q  = —A+Vq  —  A~ Aq  (2) 

where  D  —  J / At{I  +  >i+  —  A~)  with  J  a  measure  of  computational  cell  volume  and  A* 
are  the  eigenvector  split  pieces  of  the  Jacobian  matrix  according  to 

A±  =  (MT  I±T~l  M~l)A  =  A±A  (3) 

with  J*  =  ^(/  ±  sgnA)  .  Thus  as  shown  in  reference  1, 

A±Aq  =  A±AF  =  AF±  .  (4) 
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Flux  Divergence  Upwind  Methods 

To  begin  the  present  developments,  we  now  write  the  equivalent 

D  6q  =  -VF+  -  AF~  =  -A+VF  -  A~&F.  (5) 

In  two  or  more  space  dimensions  flux  difference  terms  such  as  appear  on  the  right  hand  side 
of  equation  (3)  for  each  coordinate  direction  can  be  written  either  in  strong  or  chain  rule 
(Hindman6)  conservation  law  forms  according  as  the  metric  coefficients  of  the  differential 
transformation 

dF  =  £xdF+£ydG  (6) 

are  represented  under  or  outside  of  the  difference  approximations.  In  the  latter  case,  aver¬ 
aging  the  metric  coefficients  on  the  difference  interval  is  appropriate  for  Roe  form  schemes. 
While  in  quasi  one  dimensional  flow  the  upwind  method1  with  strongly  conservative  dif¬ 
ferencing  was  found  to  be  very  robust  and  accurate,  in  2-D  or  axisymmetric  flow  the  chain 
rule  approach2  was  found  to  work  far  more  reliably.  We  will  explain  why  we  believe  this 
is  so  and  present  alternative  strong  conservation  law  formulations  that  according  to  the 
rationale  are  both  stable  and  more  accurate  than  chain  rule.  In  the  process  we  will  show 
how  to  overcome  a  problem  associated  with  nonorthogonal  grids. 

Chain  Rule  Conservation 

The  chain  rule  conservation  law  form  for  discrete  differences  is  the  natural  approxi¬ 
mation  to  equation  (6).  For  central  difference  methods  in  strong  conservation  law  form, 
it  has  long  been  recognized  and  practiced6,7'8  that  the  numerical  divergence  of  the  met¬ 
ric  coefficients  must  be  zero  under  the  operations  of  the  finite  difference  method  in  order 
for  the  method  to  meet  the  minimal  requirement  to  preserve  a  uniform  flow.  However 
in  upwind  methods,  particularly  those  in  generalized  Roe  forms,  this  fact  has  been  to  a 
greater  or  lesser  extent  ignored  and  stability  appealed  to  based  on  global  conservation  of 
the  telescoping  difference  equations.  The  problem  we  will  cite  also  exists  for  some  higher 
order  upwind  methods  with  a  numerical  flux  function  in  the  node  centered  classical  finite 
volume  cell  format  and  when  using  non  cell  face  local  metric  coefficients. 

The  issues  associated  with  local  conservation,  which  we  espouse  here,  can  be  seen 
simply  geometrically.  In  Figure  1,  we  show  the  computational  cell  associated  with  one 
sided  differencing  for  the  £  direction.  For  what  we  may  term  local  geometric  conservation 
in  the  cell,  we  consider  the  flux  vector  component  F(  for  the  £  direction  The  flux  tensor 
has  two  components  F(  and  Fn.  For  the  moment,  the  other  component  of  the  flux  tensor 
will  be  regarded  as  treated  in  another  (r))  family  of  cells  aligned  with  the  r)  coordinate 
direction.  The  flux  vector  must  be  projected  not  only  through  the  end  faces  i  and  t  +  1 
in  the  £  direction  but  also  through  the  lateral  side  faces  k  ±  The  associated  difference 
relation  is  symbolized 

A  •  •  Ff  +  6r,f)  •  (7) 


20 


1 


1 

i 

Simple  globally  conservative  one  sided  differencing  ignores  the  second  term  of  the  right  1 

hand  side  of  relation  (7)  and,  as  a  consequence,  introduces  spurious  numerical  sources  and  j 

sinks9  that  compromise  accuracy  and  stability.  1 

It  is  illuminating  and  easy  to  show  that  if  the  flux  components  for  the  lateral  r\ 
faces  of  the  cell  are  approximated  local  to  the  cell  by  |*±i=  |>  +F(  |,  +  i)a:  =  F( 

relation  (7)  is  identically  equivalent  to  chain  rule  differencing  A^F^  =  A*  F^  —  F^A^iS.  ] 

Here  we  assume  the  cell  is  properly  closed  and  the  metrics  satisfy  A(  •  S  =  0. 

Accordingly,  the  consequence  of  applying  the  chain  rule  conservation  form  is  that  in 
the  laterally  adjacent  cells  the  similar  flux  approximations  at  k  ±  1  through  the  shared 
faces  f)  ±  |  are  inconsistent  with  the  approximation  for  cell  k  and  the  method  is  in  general 
neither  locally  nor  globally  conservative.  Two  exceptions  to  the  conclusion  are  found  in 
the  cases  of  uniform  flow,  which  the  chain  rule  method  preserves,  and  no  flow  through 
the  lateral  sides  as  in  the  quasi  1-D  stream  tube  approximation  for  which  we  had  very 
good  success1  in  strongly  conservative  differencing.  Thus  in  principle  the  chain  rule  and 
simple  globally  conservative  locally  one  dimensional  methods  are  relatively  more  effective 
where  the  grid  is  nicely  flow  aligned,  smooth  gradients  are  weak  and  discontinuities  are 
orthonormal.  In  practice  chain  rule  conservation,  with  similar  one  sided  difference  terms 
for  the  t)  coordinate  direction,  works  quite  well2,3  for  the  numerically  dissipative  first  order 
method.  For  second  order  methods  with  limiters  the  error  in  conservation  is  diminished 
but  the  methods  also  decline  in  robustness,  about  which  more  will  be  said  later. 

Locally  One  Dimensional  Strong  Conservation 

The  conclusion  from  the  above  discussion  is  that  a  method  closely  akin  to  the  chain 
rule  locally  one  dimensional  one  sided  upwind  differencing  procedures  of  references  2,  3  and 
10  but  which  features  consistent  local  geometric  conservation  for  each  flux  component  in 
its  own  family  of  one  sided  cells  is  obtained  by  replacing  the  simple  one  sided  flux  difference 
A  {F{  for  equation  (5)  by  cell  flux  divergences  symbolized  by  A t  •  F$  with  the  lateral  fj 
side  fluxes  in  equation  (7)  given  a  symmetric  definition,  for  example  averaged  between  k 
and  k  +  1.  Similar  generalized  one  sided  upwind  difference  expressions  are  given  for  the  77 
coordinate  direction  and  we  can  write 

DSq  =  -A+V*  'Ft-AJAfFt-  A+V„  .  F,  -  A+An  .  (8) 

This  approach  is  both  globally  conservative  and  more  consistent  with  the  accuracy 
aims  of  local  conservation  and  accordingly  is  expected  to  be  more  stable  and  robust. 

The  above  locally  one  dimensional  approach  has  a  numerical  virtue  in  that  the  splitting 
operators  apply  homogeneously  and  rationally  to  their  own  convective  fluxes. 

There  are,  however,  drawbacks  still  to  be  addressed.  The  first  of  these  deals  with 
nonorthogonal  meshes.  While  the  problem  is  not  observed  numerically  to  be  severe  for 
moderate  deviations  from  orthogonality,  from  equation  (8)  it  is  obvious  that  in  the  limit 
as  the  mesh  is  skewed  and  the  rj  and  £  fluxes  become  identical,  other  effects  aside,  the 
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method  loses  information  from  one  sector  of  the  flow  domain  and  its  solution  must  be¬ 
come  inconsistent  with  the  multi-dimensional  PDE.  A  second  weakness  is  that  for  general 
curvilinear  coordinates  the  £  direction  in  the  Cartesian  bases  changes  from  cell  to  cell  so 
there  is  not  a  unique  and  consistent  direction  on  which  to  effect  the  flux  decomposition 
and  satisfy  assumptions.  A  third  weakness  is  that  in  the  locally  one  dimensional  method 
neither  the  partial  nor  total  flux  is  identically  conserved  locally  in  each  family  of  cells.  The 
significance  is  that  residual  reduction  at  convergence  may  not  be  as  firm  as  with  the  last 
method. 

Symmetric  Locally  One  Dimensional  Conservation 

Objections  one  and  two  above  may  be  overcome  by  constructing  methods  that  are 
symmetric  overall  in  the  two  families  of  cells.  By  this  we  mean  that  all  the  cells  of  the 
procedure  are  shared  among  the  two  coordinate  directions;  a  £  one  sided  difference  cell  for 
one  coordinate  is  also  an  77  cell  for  the  other  coordinate.  A  last  requirement  for  symmetry 
is  that  a  unique  definition  of  total  flux  be  given  for  every  cell  face.  These  conditions  ensure 
that  regardless  of  the  bases  for  flux  tensor  decomposition  in  adajacent  cells,  the  total  flux 
that  numerically  flows  out  of  a  face  of  one  cell  also  flows  through  the  same  face  into  its 
neighbor  cell. 

Beyond  the  classical  node  centered  finite  volume,  two  kinds  of  computational  cells 
can  lend  themselves  to  the  construction  of  symmetric  strongly  conservative  biased  upwind 
methods.  One  of  these  cell  types  is  the  natural  mesh  cell  bounded  at  the  four  corners 
(in  2-D)  by  its  data  nodes.  The  other  is  the  laterally  symmetric  and  upwind  (coordinate) 
biased  cell  of  Figure  1. 

Symmetric  Methods  in  Natural  Mesh  Cells 

Of  two  approaches  we  have  considered,  one  splits  the  weighted  average  of  the  flux 
divergences  in  adajacent  mesh  cells  on  either  side  of  the  differencing  coordinate  line  with 
the  eigen  state  determined  by  adjacent  node  averages  on  the  difference  interval.  This 
scheme  has  laterally  wide  support  in  the  difference  stencil  and  thus  is  not  very  accurate  in 
high  gradient  regions  and  will  smear  solutions.  The  other  approach  splits  the  side  averaged 
flux  difference  according  to  the  eigen  state  at  the  mesh  cell  centers.  This  procedure  implies 
compact  differencing  schemes  and  is  restrictive.  The  method  also  requires  a  lot  of  averaging 
which  can  be  mitigated  by  placing  data  nodes  at  the  cell  centers.  The  latter  idea  is  not 
bad  as  we  shall  consider  next. 

Symmetric  Methods  on  Biased  Upwind  Cells 

Additional  data  nodes  at  the  (body)  centers  of  the  natural  mesh  cells  create  a  sym¬ 
metric  geometric  pattern  in  which  the  biased  upwind  cells  (Figure  2)  of  the  base  level  mesh 
for  one  coordinate  direction  are  the  similarly  biased  upwind  cells  in  the  other  coordinate 
direction  of  the  body  centered  mesh.  Since  the  nodes  of  the  body  centered  mesh  are  cen¬ 
tered  on  the  lateral  faces  of  the  base  mesh  upwind  cells  and  visa  versa,  there  is  a  unique 
definition  of  the  total  flux  at  the  node  on  every  cell  face. 


For  each  face  of  every  cell  the  flux  tensor  can  be  decomposed  in  orthogonal  components 
according  to  reference  directions  that  may  be  determined  independently  from  cell  to  cell. 
The  orthogonal  decomposition  can  be  expressed  in  the  unit  base  vectors  a  and  0  as 

F  =  aa»F  +  00»F  =  Fa+F^  (9) 

We  assume  the  base  vectors  are  related  to  and  generally  align  with  the  coordinate  lines. 
Then  for  either  component  we  can  construct  the  locally  one  dimensional  flux  divergence 
Aa  •  Fa  according  to  equation  (7)  with  Fa  replacing  F (.  Of  course  a  similar  relation  is 
developed  for  A p*F^  appropriately  interchanging  the  upwind  one  sided  and  lateral  central 
difference  symbols  among  £  and  q. 

For  either  flux  divergence  component  and  by  analogy  with  the  chain  rule  CSCM 
method,  generalized  flux  difference  eigenvector  splitting  may  be  carried  out  based  on  the 
associated  upwind  difference  interval  averaged  eigenstate  to  give 

A  •  F±  —  v4(/±)A  •  F  (10) 

Finally  associating  coordinate  £  (j)  with  the  a  flux  tensor  component  and  r?  (k)  with  0, 
we  can  write  the  flux  divergence  split  implicit  method  as 

[A{ii  +  a+)v<  +  A(J-  +  a;)a<  +  B(i;  +  A+)V,  +  b(i;  +  a;)a„]%,* 

=  -(i(7+)AJ_1  •  F°  +  A(f~) Ay  •  Fa  +  5(/;)A,_1  .  F0  +  B(I~)Ak  .  Fp)  (11) 

Here,  for  a  two  level  linearized  method,  we  have  via  Eq.  (5)  approximated  the  split 
Jacobian  matrix  A±  by  the  similarity  transform  based  on  the  A  operator  and  the  associated 
difference  interval  estimated  eigenvalues 

A*  cs  [W*  +  lc  (12) 

For  very  skewed  meshes,  a  better  estimate  of  the  eigenvalues  would  be  required. 

While  inherently  more  accurate,  the  symmetric  staggard  grid  method  described  above 
is  seen  to  closely  approximate  the  chain  rule  CSCM  algorithm.  We  have  given  the  new 
derived  method  the  name  S-CSCM. 

Local  Total  Flux  Strong  Conservation 

Finally  we  give  a  simple  approach  which  provides  total  flux  conservation  in  each 
computational  cell.  In  the  approach  we  replace  in  Eq.  (5)  the  partial  fluxes  by  the  total 
flux  to  give 


D6q  =  .  F  -  A- Ae  .  F  -  /i  +  V„  .  F  -  A~  An  .  F  (13) 

With  Eq.  (13)  a  sufficient  (and  we  would  guess  necessary)  condition  for  convergence 
to  steady  state  is  that  the  flux  divergence  vanish  for  every  computational  cell.  Within 


23 


the  total  flux  divergence  framework  the  subscripts  on  the  symbolic  divergence  operators 
Vj  •  F  only  serve  to  identify  a  kind  of  cell  in  the  grid.  Further,  the  subscripts  £  and  r? 
on  the  projection  operators  serve  to  identify  characteristic  directions  in  which  the  velocity 
vectors  will  be  composed  in  normal  and  tangential  components  for  advection.  When  these 
directions,  which  we  can  free  from  constraint  to  either  or  both  the  curvilinear  coordinate 
directions,  are  taken  to  be  orthogonal,  it  can  be  shown  that  within  truncation  errors  similar 
to  those  for  locally  one  dimensional  conservative  eigenvector  split  upwind  differencing  the 
method  is  consistent  with  the  PDE. 

Based  in  an  older  and  independently  conceived  notion  by  Roe  and  Lombard  the  local 
cell  total  flux  divergence  A  •  F  can  be  identified  with  a  residual  (“fluctuation”,  Roe)  that 
is  “sent”  along  grid  lines  by  the  method  through  the  projection  operators  Af>r}  to  relax  the 
state  at  appropriate  upwind  mesh  points.  With  freedom  of  the  characteristic  directions 
from  the  grid,  we  achieve  greater  freedom  in  how  and  where  to  send  residual  information. 

Second  Order  TVD  Schemes 

For  the  new  locally  one  dimensional  procedure,  higher  order  spatial  difference  meth¬ 
ods  for  the  right  hand  side  of  equation  (13)  can  be  constructed  of  the  split  partial  flux 
divergence  by  pieced  linear  combinations  of  the  standard  upwind  and  central  difference 
correction  operators.  In  reference  10,  Lombard  proposed  a  simple  limiter  (Scheme  Illb) 
for  the  conservative  flux  difference  split  pieces  of  such  a  second  order  method.  That  lim¬ 
iter  by  direct  analogy  can  be  used  here  to  enhance  nonlinear  stability.  We  note  that  the 
natural  second  order  method  for  a  scheme  constructed  of  the  total  flux  divergence  involves 
the  “Fromm”  weights  of  half  central  and  half  upwind. 

Viscous  Terms 

For  both  the  locally  one  dimensional  and  total  flux  formulations  it  is  feasible  following 
the  work  of  reference  2  to  simply  add  viscous  terms  in  a  conventional  thin  layer  central 
difference  formulation.  For  the  locally  one  dimensional  approach  this  may  be  appropriate. 
But  from  the  strong  conservation  point  of  view  such  a  practice  will  interfere  with  the 
desirable  feature  of  zeroing  the  total  flux  divergence  in  each  computational  cell-  To  preserve 
that  property  of  the  inviscid  method,  we  have  only  to  add  the  viscous  stress  terms  at  the 
boundaries  of  the  cells  such  that  the  terms  contribute  consistently  to  the  total  flux  whether 
or  not  it  is  subsequently  decomposed  into  locally  one  dimensional  components. 

Numerical  Results 

We  show  here  results  of  computations  for  a  topical  difficult  propulsion  base  flow  prob¬ 
lem  of  a  wind  tunnel  model  tactical  missile  flow  involving  realistically  complex  geometry 
and  flow  structure  including  a  centered  propulsive  jet.  The  freestream  Mach  number  is 
1.4,  base  diameter  Re  0.92x10°  and  nozzle  exit  Mach  number  2.7  and  static  pressure  ratio 
of  2.5.  The  problem  has  recently  been  solved  with  a  geometry  flexible  multiple  patched 
mesh  CSCM  code  -  based  on  the  chain  rule  conservation  formulation. 
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The  problem  was  solved  previously  by  us  with  an  adapted  single  block  mesh  code 
on  a  double  wraparound  grid  of  Figure  3.  The  grid  generated  in  patches  with  a  simpli¬ 
fied  algebraic  procedure  has  nice  quality  over  the  external  and  internal  corners.  But  the 
constrained  topology  is  seen  to  waste  mesh  points  through  excessive  clustering  in  the  wake. 

The  new,  better  balanced  mesh  for  the  present  solutions  is  shown  in  Figure  4.  A  plot 
of  the  outline  of  the  multiple  patch  structure  of  the  mesh  in  the  vicinity  of  the  nozzle  exit 
is  shown  in  Figure  5.  In  the  boundary  layer  region  the  mesh  has  an  effective  wraparound 
character  on  an  overset  patch  local  to  the  corner. 

We  show  in  Figure  6  a  Mach  contour  plot  for  the  solution  on  the  new  mesh.  On  the 
whole  the  better  balanced  mesh  employing  about  the  same  number  of  points  produces 
a  result  of  crisper  flow  structural  detail  than  the  less  flexible  topology  whose  solution  is 
shown  in  Figure  7.  Nicely  evident  in  Figure  6  is  a  weak  lip  shock,  shear  layer,  barrel  shock, 
(near)  contact  discontinuity  and,  above  them  all,  a  recompression  shock  in  the  external 
flow. 

Lastly  in  Figure  8,  we  show  a  comparison  of  computational  results  and  experiment. 
The  wave  in  the  base  region  is  due  to  two  counterrotating  persistent  vortices  in  the  base 
bubble.  We  note  at  the  corners  a  strong  qualitative  disagreement  between  solutions  on 
the  wraparound  grids  and  a  third  nonwraparound  step  mesh13.  Evidently  there  the  mesh 
topology  strongly  influences  the  mechanism  of  separation.  This  point  requires  further 
study  into  topology  and  mesh  refinement. 
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5.  Uniformly  Second  Order  Accurate  ENO  Schemes  for  the 

Euler  Equations  of  Gas  Dynamics 

Very  recently,  a  new  class  of  uniformly  high  order  accurate  essentially  non-oscillatory 
(ENO)  schemes  have  been  developed  by  Harten  and  Osher1  and  Harten  et.  al.2~4.  They 
presented  an  hierarchy  of  uniformly  high  order  accurate  schemes  which  generalize  Go¬ 
dunov’s  scheme5,  and  its  second  order  accurate  MUSCL  extension6,7  and  TVD  scheme8 
to  arbitrary  order  of  accuracy. 

In  contrast  to  the  earlier  second-order  TVD  schemes  which  drop  to  first-order  accu¬ 
racy  at  local  extremums  and  maintain  second-order  accuracy  at  smooth  regions,  the  new 
ENO  schemes  are  uniformly  high  order  accurate  thoughout  even  at  ciitical  points.  The 
ENO  schemes  use  a  reconstruction  algorithm  which  is  derived  from  a  new  interpolation 
technique  that  when  applied  to  piecewise  smooth  data  gives  high  order  accuracy  whenever 
the  function  is  smooth  but  avoids  a  Gibbs  phenomenon  at  discontinuities.  An  adaptive 
stencil  of  grid  points  is  used;  therefore,  the  resulting  schemes  are  highly  nonlinear  even  in 
the  scalar  case. 

Theoretical  results  for  the  scalar  constant  coefficient  case  and  numerical  results  for 
the  scalar  conservation  law  and  for  the  one-dimensional  Euler  equations  of  gas  dynamics 
have  been  reported  with  highly  accurate  results.  Such  high  order  ENO  schemes  have  the 
potential  to  be  adapted  to  the  current  Euler/Navier-Stokes  flow  solvers  as  one  does  for  the 
TVD  schemes  to  further  enhance  the  accuracy  of  flowfield  simulation.  Implementation  can 
be  either  as  a  genuine  higher-order  flow  solver  as  in  the  present  work  or  as  a  post-processor 
to  enhance  the  resolution. 

In  this  paper  we  first  formally  extend  our  second-order  TVD  schemes  described  in9-13 
to  uniformly  second-order  ENO  schemes  for  the  two-dimensional  Euler  equations  in  the 
general  coordinate  systems.  Both  explicit  and  implicit  schemes  are  described.  It  is  shown 
that  our  previous  TVD  schemes  are  a  special  case  of  the  present  ENO  scheme.  Here,  the 
TVD  requirement  is  replaced  by  a  less  restricted  essentially  non-oscillatory  condition,  a 
concept  advanced  by  Harten  et  al1-4. 

We  apply  the  resulting  schemes  to  simulate  complex  2-D  gas  flows  involving  multiple 
discontinuities  to  test  the  global  accuracy  of  ENO  schemes.  In  particular,  we  compare  the 
ENO  results  directly  with  the  TVD  results  for  the  same  simulation;  thus  merits  of  each 
method  can  be  assessed. 

In  Section  2,  the  two-dimensional  Euler  equations  in  the  general  coordinate  system 
are  formulated  in  the  framework13  of  a  characteristic  flux  difference  split  method.  Such 
methods  that  we  employ  are  CFDS14  and  CSCM15.  Properties  of  the  Euler  equation 
which  are  relevant  to  the  present  development  are  briefly  reviewed.  To  achieve  higher 
order  accuracy,  a  modified  flux  similar  to  Harten ’s  is  introduced. 

In  Section  3,  we  describe  a  modified  flux  which  can  have  either  an  ENO  or  a  TVD  prop¬ 
erty.  The  relationship  between  the  present  approach  and  Harten-Osher’s  UNO  scheme  is 
illustrated.  Both  explicit  and  implicit  schemes  are  described  for  the  2-D  Euler  equations  in 
general  coordinate  systems.  Strang-type  dimensional  splitting16  is  used  for  explicit  meth¬ 
ods  to  simulate  time-dependent  problems.  Standard  ADI17,18  and  diagonally  dominant 
ADI13,15  solution  procedures  for  steady-state  calculations  are  discussed. 
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Numerical  experiments  with  present  ENO  schemes  for  a  1-D  shock  tube  problem, 
truly  nonstationary  shock  wave  diffraction  around  a  cylinder,  shock  wave  collision  over  a 
circular  arc,  and  a  steady  transonic  flow  over  a  circular  arc  in  a  channel  are  reported  in 
Section  4.  Numerical  results  are  compared  with  either  analytic  or  experimental  results  to 
assess  the  present  work. 

Conclusions  and  suggestions  are  given  in  Section  5. 

Theoretical  Considerations 

The  governing  equations  of  two-dimensional  unsteady  gas  dynamics,  neglecting  the 
effects  of  viscosity  and  heat  transfer,  in  general  curvilinear  coordinate  systems  (£,77)  arc 
written  in  conservation  form  as 

dTQ  +  dzF  +  dnG  =  0  (1) 

where  Q  =  Q/J  and  F  =  (£<Q  +  £XF  +  iyG)/J,  G  -  ( r)tQ  +  rjxF  +  r]yG)/J,  and  J  = 
tx Vy  —  £yrix,  the  metric  Jacobian.  Q  —  (p,pu,pu,  E)r  is  the  conservative  variables,  F  = 
(pu, pu2  +  p, puv,u[E  +  p))T  and  G  =  ( pu, ptm , pu2  +  p,  v(E  +  p))T  are  the  flux  vectors  in 
Cartesian  coordinates. 

Here  p  is  the  gas  density,  u,  v  are  the  gas  velocity  components,  E  the  total  energy 
and  p  the  gas  pressure.  The  pressure  is  related  to  other  variables  by  the  equation  of  state, 
for  a  perfect  gas,  p  =  (7  —  1)[.E  —  0.5p(u2  +  v2)],  where  7  is  the  specific  heats  ratio. 

Due  to  the  hyperbolicity  of  system  (l),  the  Jacobian  coefficient  matrix  ^  =  dF/dQ 
of  the  transformed  equations  has  real  eigenvalues 

=  U,  a2  =  U  +  c^,  <13  =  U,  and  a4  =  U  —  (2) 

with  U  =  &  +  and  =  c^J £2  -f-  £2,  where  c  =  ^7 p/p,  is  the  speed  of  sound. 

One  can  also  find  similarity  transformation  matrices  T(  and  Tn  that  diagonalize  the 
system 

T;lA&  =  At  =  diag{a;}, 

T~1Br]T,j  =  A.„  =  diag{6(}  (3) 

Since  Eq.(l)  is  a  hyperbolic  conservation  law,  both  features  of  hyperbolicity  and  con¬ 
servation  can  be  utilized  in  constructing  numerical  methods  for  solving  them.  A  simple 
and  natural  way  to  unify  these  two  aspects  is  to  put  Eq.(l)  in  the  following  form: 

drQ  +  (A+  +  AJ)dzF  +  (B+  +  B-)dnG  =  0  (4) 

where  and  B *  are  the  split  normalized  Jacobian  coefficient  matrices  that  effect  the 

flux  difference  splitting13.  Analytic  expressions  for  A *  and  B *  are  given  in  Appendix  1. 

Several  ways  of  constructing  higher-order  schemes  based  on  (4)  have  been 
investigated9-11. 
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One  way  to  construct  a  family  of  second  order  accurate  schemes  is  to  approximate 
Eq.(4)  by 

drQ  +  (i+  +  AJ)d(FM  +  (JB+  +  B~)dnGM  =  0  (5) 

where  FKt  =  F  +  E  and  GM  =  G  +  H  are  modified  fluxes. 

In  the  following  based  on  Eq.(5),  a  numerical  method  of  uniformly  second  order  ac¬ 
curacy  in  time  and  space  which  combines  both  characteristic  and  conservation  features  of 
Eq.(l)  is  described. 

Uniformly  Second  Order  Accurate  ENO  Schemes 

Let  us  define  a  uniform  computational  mesh  system  (fy,??*)  with  mesh  sizes  A£,  and 
Atj  and  let  Q™ k  denote  the  value  of  Q  at  time  level  nAr  and  at  position  (jA^,kArj). 

A  conservative  difference  scheme  for  Eq.(5)  can  be  expressed  in  terms  of  numerical 
fluxes  Fn  and  GN  as  follows: 


(6) 


In  the  following,  only  the  L(  operator  is  given  in  detail.  Similar  expression  can  be 
given  for  the  Ln  operator.  Explicit  schemes  for  Eq.(6)  can  be  constructed  by  time  splitting. 
For  the  L ^  operator  in  the  ^-direction,  we  have 


£<<3?.*  =  Qh  =  Q?.k  -  McAi,*  - 

where  A ^  =  Ar/A£ 

In  Eq.(7),  FV+^  fc  is  the  numerical  flux  which  is  given  by 


,  =  fM_,  k  —  At  .  A,**  ,Fm 


]  +  \,k  r]  +  l,k  A£;  +  i+|.* 


with  Ay+|()  =  ()y+x  -  Oy. 

In  Eq.(8),  the  split  “normalized”  Jacobian  coefficient  matrix  A ^  satisfies 


(7) 


(8a) 

(86) 


T-'AfTe  =  If  =  diag  {a?}, 


1  ±  sgn(a;) 


(9) 


The  value  of  E  at  nodal  point  j,k  for  the  present  ENO  scheme  is  E:-,k  = 
(e i ,  e2, ...,  64)  Jk  and  its  l  components 


e; 


-pm(A+elj+ik,A-er+ik),el}_ik 
+  (3m(A+el]_i  h,A-elj_i  J] 


(10) 


Here  m  is  the  min  mod  function: 

m(a,  b)  =  smin(|a],  |6|)  if  sgn  a  =  sgn  b  =  s 

=  0,  otherwise  (12) 

The  m  function  is  defined  as: 

m(a,  b)  =  a  if  |a|  <  |6| 

=  b  if  |a|  >  |6|  (13) 

For  0  ~  0.,  we  have  a  second-order  TVD  scheme  like  before.  For  0  =  0.5,  we  have  an 
uniformly  strictly  non-oscillatory  scheme. 

For  scalar  constant  coefficient  case,  e.g.  dtu  +  adxu  =  0,  Eq.(10)  can  be  proven  to  be 
identical  to  Harten-Osher’s  non-oscillatory  MUSCL-type  scheme1  which  uses  reconstruc¬ 
tion  via  deconvolution  (RD)  approach,  here  of  degree  2 
That  is 

R(x\vn)  —vf  +  (x  —  xy)A  j+^v/h 

+  \(x~  xj)(x  ~  *y+i)(«/+i  -  sj)/h 
for  Xj  <  x  <  xj+ 1  (14) 

with 

Sj  =  ~  0,vn^'7iH^Xj' +  0;  (15) 

where  /^(xjv)  is  a  non-oscillatory  piecewise  parabolic  interpolation. 

To  further  illustrate  the  present  approach  and  its  relationship  to  Harten-Osher’s  UNO 
scheme1,  let  us  consider 

dtu  +  adxu  =  0,  a  >  0 
Then  from  Eqs.(7),  (8)  and  (10)  we  have 

v;+1  = »;  -  7  - 

+  iBS[A.(i^)A-r?  ,  A.(i^)A+v»]  .  (l^)A  +  v? 

-  im(A.(i^)A+v?  ,  A+(i^)A  +  v“]})  (16) 
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Here  a  =  \a  and  it  can  be  a  variable.  For  a  =  1,  a  =  A,  this  reduces  to  Eq.(4.1G)  of 
reference  3. 

It  is  rather  constructive  to  have  the  scheme  written  in  this  form,  since  one  can  easily 
extend  that  for  a  single  wave  equation  to  a  scalar  conservation  law  and  to  hyperbolic 
systems  of  conservation  laws  in  a  straightforward  manner. 

One  can  similarly  construct  schemes  of  higher  than  second  order  accuracy  for  the 
hyperbolic  system  of  conservation  laws.  For  example,  in  reference  3,  the  scheme  defined 
by  Eq.(4.19)  was  obtained  using  N  =  3  in  the  HP  approach. 

In  Eq.(lO),  one  can  replace  both  of  the  m  by  m  to  obtain  another  truly  non-oscillatory 
scheme1.  For  further  details  of  the  method,  the  reader  is  encouraged  to  read  the  original 
papers1-4. 

It  is  noted  that  the  present  approach  uses  a  high  order  representation  of  the  flux  vector, 
in  contrast  to  the  approach  defined  by  Eq.(14)  which  uses  a  high  order  interpolation 
for  the  conservative  variable  Vj.  As  a  result,  the  present  approach  is  quite  different  from 
that  of  Harten  et  al1-4  for  nonlinear  cases. 

Implicit  Schemes 

Implicit  schemes  can  be  consistently  constructed  based  on  the  explicit  schemes.  An 
implicit  scheme  using  back  Euler  in  time  can  be  given  as 

[/  +  i+A_AM  +  A~A  +  Am  +  £+A_£M 

+  B~A  +  BKf}6Qn  =  RHS  of  Eq.  (6)  (17) 

where  AM  -  dFM  jdQ  and  BM  =  dGM /dQ. 

In  practice  we  adopt  the  following  approximation 

Am  =  A  +  dE/dQ  «  A 

and  Bm  =  B  +  dH/dQ  «  D  (18) 

This  corresponds  to  a  linearization  procedure  and  reduces  the  computational  efforts  con¬ 
siderably. 

Block  Tridigonal  Approximate  Factorization 

The  LHS  of  equation  (17)  may  be  approximately  factored  in  two  forms 
Standard  ADI  (Beam  and  Warming17,18) 

\-A+,D(,A-}\-B+,Dn,B-}6Q  =  RHS  (19a) 

where  D (  —  I  - 1-  |A|,  Dn  =  I  +  |£|  or  with  reduced  errors  by 
DDADI  (Lombard  et  al15) 

\-A  +  ,D,A-}D-l\-B+,D,B-}6Q  (196) 

with  D  =  I+\A\  +  \B\. 
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Equations  (19)  gives  the  block  tridigonal  solution  sequences 

=  RHS 

[-B+,Dn,B-]£Qn  =  6Q * 

Qn+1  =  Qn  +  6Qn  (20a) 

and 

[-A+,D,A~\6Q*  =  RHS 
l-B+,D,B~]6Qn  =  D6Q * 

Qn+1  =  Qn  +  6Qn  (20  b) 

Implicit  TVD  schemes  with  DDADI  and  bidiagonal  approximate  factorization  have  been 
reported  by  Yang  et  al.13. 

Numerical  Results  and  Discussion 

Results  of  several  numerical  experiments  are  given  in  the  following.  Both  unsteady 
flow  using  the  explicit  method  and  a  steady-state  solution  using  the  implicit  method  are 
included. 

(1)  1-D  shock  tube  solution 

Fig.  1  shows  the  solution  of  the  1-D  shock  tube  problem  of  Sod19.  Here  we  have  used 
Eq.(lO)  with  both  m  replaced  by  minmod  function  m.  The  same  results  using  a  TVD 
scheme  are  shown  in  Fig.  2  for  comparison.  Moderate  improvement  at  the  expansion  fan 
and  the  contact  surface  is  observed.  Slight  attenuation  is  observed  for  the  ENO  result 
which  is  allowed  by  the  oscillation  diminishing  property. 

(2)  Shock  wave  reflection  by  a  circular  cylinder 

Fig.  3  displays  the  density  contours  and  Mach  number  contours  for  a  Mach  shock 
number  M,  =  2.81  using  the  ENO  scheme.  A  comparison  with  an  experimental  Schlieren 
picture  is  shown  in  Fig.  4.  Excellent  agreement  is  found  in  nearly  every  aspect,  except  for 
some  viscous  effects  not  accounted  for  by  the  Euler  equations. 

(3)  Shock  wave  collision  over  a  15%  circular  arc 

Fig.  5  shows  the  density,  pressure  and  Mach  number  contours  for  two  shock  waves 
(with  M t  =  10,  7  =  1.1)  in  a  head  on  collision  over  a  circular  arc  bump  using  the  ENO 
scheme.  The  flowfield  is  very  symmetric.  The  sudden  jump  in  the  pressure  after  the 
collision  can  be  clearly  seen.  Multiple  Mach  shocks  (triple  shock)  and  sliplines  and  their 
interaction  such  (vortices)  are  accurately  resolved.  Attached  leading  edge  bow  shock  is 
observed  in  contrast  to  a  detached  bow  shock  for  7  =  1.4. 

(4)  Asymmetric  Shock  Wave  Collision  (7  =  1.4) 

Figure  6  displays  the  post  collision  flow  of  a  shock  wave  moving  from  the  right  to 
the  left  with  a  delayed  shock  moving  from  the  left  to  the  right.  A  vertical  contact  surface 
behind  the  left  moving  shock  (about  to  exit  the  frame)  and  curved  bow  shocks  reflected 
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from  either  side  are  prominent  features  of  the  solution.  Complicated  flow  structures  such 
as  vortices,  multiple  Mach  shocks  and  sliplines  can  be  identified. 

The  above  numerical  experiments  demonstrate  the  capability  of  the  present  ENO 
scheme  for  time-dependent  problems;  the  quality  of  results  is  quite  good. 

Next  we  consider  a  steady-state  aerodynamic  problem  using  an  implicit  ENO  scheme. 

(5)  GAMM  transonic  channel  flow 

We  considered  one  of  the  standard  test  cases  given  at  the  GAMM  Workshop,  namely, 
an  internal  two-dimensional  transonic  flow  through  a  parallel  channel  having  a  4.2%  thick 
circular  arc  at  the  lower  wall.  The  ratio  of  static  downstream  pressure  to  total  pressure 
is  0.623512,  corresponding  to  =  0.85  in  the  isentropic  flow.  The  Cp  distribution  at 
the  lower  surface  and  the  Mach  number  contour  plot  for  the  DDADI  solution  with  ENO 
and  TVD  schemes  are  shown  in  Fig.  7  and  Fig.  8,  respectively.  Both  are  calculated  with 
CFL=10.  The  ENO  scheme  exhibits  superior  results  to  those  of  the  TVD  scheme. 

Concluding  Remarks 

Uniformly  second-order  accurate  explicit  and  implicit  schemes  with  the  essentially 
non-oscillatory  property  have  been  described  for  the  two-dimensional  Euler  equations  in 
general  curvilinear  coordinates.  Applications  to  time-dependent  and  asymtotic  steady- 
state  aerodynamic  problems  have  been  carried  out.  The  results  exhibit  robust  stability 
and  high  resolution  of  flowfields  involving  multiple  discontinuities.  In  comparison  with 
second-order  TVD  schemes,  slight  to  considerable  gain  in  performance  can  be  achieved 
for  a  rather  small  increase  in  computational  effort  depending  upon  the  complexity  of  the 
flowfields  involved.  One  of  the  main  objectives  is  to  test  the  new  ENO  scheme,  which  has 
a  completely  different  underlying  design  principle  in  practical  aerodynamic  applications. 
Numerical  experiments  suggest  that  potential  benefit  can  be  gained  with  this  new  class 
of  high  order  essentially  non-oscillatory  schemes.  Extension  to  include  the  viscous  terms 
for  the  Navier-Stokes  equations  can  be  accomplished  in  a  standard  manner  using  central 
differencing.  Such  implementation  is  currently  underway. 
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Analytical  expressions  for  A±  and  B 
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The  expressions  for  5*  are  obtained  by  replacing  £x 
and  £y  by  fjx  and  fjy,  respectively. 
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Figure  5.  Time  dependent  symmetric  shock  diffraction  - 
collision  over  circular  arc  airfoil,  (a)  Density  contours, 
(b)  Pressure  contours,  (c)  Mach  number  contours. 
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Figure  6.  rime  dependent  asymmetric  shock  diffraction  -  collision 
over  c i rc>  iar  arc  airfoil,  (a)  Density  contours,  (b)  Pressure 
contours,  (c)  Mach  number  contours. 
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6.  Three  Dimensional  Hypersonic  Shock  Capturing  CSCM  Upwind 
Implicit  Navier-Stokes  Method 

The  effective  engineering  design  of  the  new  transatmospheric  space  vehicles  flying  at 
hypersonic  speeds  requires  computational  aerothermodynamics  (CAT)  technology  in  the 
frontier  of  present  scientific  research1.  The  aeroassisted  orbital  transfer  vehicle  (AOTV)  to 
fly  at  Mach  30  in  the  upper  atmosphere  has  the  design  requirement  to  predict  the  large  low 
speed  zone  of  recirculating  flow  with  the  turbulent  wake  behind  this  vehicle  and  represents 
a  major  research  activity.  Another  major  hypersonic  vehicle  under  current  investigation 
is  the  transatmospheric  vehicle  (TAV)  which  will  fly  into  orbit  from  Earth’s  surface  and 
land  back  using  its  own  power.  The  design  of  these  vehicles  will  require  the  integration 
of  different  areas  of  technology.  In  particular,  the  combination  of  an  efficient  three  di¬ 
mensional  shock  capturing  CFD  method  such  as  the  3-D  CSCM  upwind  method2  with 
an  efficient  nonequilibrium  chemistry  procedure3  is  needed  to  account  for  the  significant 
effects  of  thinner  shock  layer  on  drag  and  wall  heat  transfer  associated  with  high  altitude 
hypervelocity  flight. 

The  3-D  CSCM  algorithm  employed  here  is  an  implicit  method  of  planes  symmetric 
Gauss-Seidel  relaxation  scheme4.  The  data  is  conveniently  stored  on  successive  planes 
along  the  streamwise  coordinate,  and  the  system  of  equations  is  iterated,  i.e.  solved  twice, 
in  marching  successive  planes  of  the  streamwise  coordinate;  along  the  forward  direction, 
first,  and  along  the  backward  direction,  afterwards.  In  each  plane  the  solution  is  ob¬ 
tained  by  a  two  level  pseudo  time  dependent  relaxation  procedure  based  on  the  diagonally 
dominant  approximate  factorization  DDADI5.  The  space  marching  alternating  directional 
sweeps  in  the  streamwise  coordinate  are  von  Neumann  unconditionally  stable6,7,8  for  zones 
of  subsonic  and  streamwise  separated  and  reversed  flows  as  well  as  supersonic  flow.  Much 
as  do  the  more  restrictive  PNS  techniques,  the  present  space  marching  method  results  in 
improved  propagation  of  nonlinear  advection  to  accelerate  convergence  to  steady  state. 

The  relaxation  algorithm  has  half  the  storage  requirements  and  substantially  acceler¬ 
ated  convergence  over  related  two  level  linearized  time  dependent  implicit  methods.  The 
method  also  requires  only  a  few  2-D  cross-flow  data  planes  in  core  at  any  marching  step 
and  thus  with  fast  disc  storage,  such  as  is  available  on  the  CRAY-XMP  machines,  can 
efficiently  treat  complex  problems  requiring  very  large  numbers  of  mesh  points  in  3-D. 
Accordingly  the  method  provides  a  base  where  chemistry  procedures  can  efficiently  be 
exploited. 

The  purpose  of  this  present  report  is  to  describe  important  features  of  the  3-D  CSCM 
shock  capturing  method  for  hypersonic  flow,  and  to  show  the  performance  of  the  method 
in  the  simulation  of  the  hypersonic  equilibrium  flowfield  around  a  sphere-cone-cone-flare 
reentry  vehicle  (RV)  at  angle  of  attack.  The  present  work  for  perfect  gases  has  recently 
been  extended  to  flows  of  general  equilibrium  chemically  reacting  gases  by  Nagaraj,  et  al9. 
The  3-D  code  is  substantially  validated  by  comparison  with  results  in  axisymmetric  flow 
for  the  generically  similar  2-D/axisymmetric  CSCM  upwind  code  that  has  previously  been 
extensively  tested  against  experiment  in  both  external  and  internal  flows  as  described  in 
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the  referenced  works.  For  both  classes  of  flows  CSCM  has  demonstrated  convective  heat 
transfer  capability10. 

In  the  following  sections  we  will  first  describe  the  numerical  method  and  then  present 
the  numerical  results. 

CSCM  Gasdynamic  Formulation 

The  conservative  finite  volume  finite  difference  formulation  of  the  3-D  compress¬ 
ible  Navier-Stokes  equations  in  general  curvilinear  coordinates  i[x,y,z),  r](x,y,z),  and 
0(x,  y,  z)  is  expressed  as 

Jdtq+l&F  +  JydtG  +  t,dtH+ 

T)xdr>F  +  VydnG  +  rjzdnH+ 

rpxd^,F  +  rpyd^G  +  tpzd^H  =  0  (l) 

where  F,  G,  and  H  are  the  flux  vectors  in  the  Cartesian  coordinate  directions  x,  y,  z, 
respectively;  while  q  represents  the  conservative  dependent  variables, 

q  =  (p  ,  pu  ,  pv  ,  pw  ,  e)  (2) 

F  =  [pu  i  pu2  +p-  axx  ,  puv  -  Txy  ,  pu w  -  txx  , 

(e  +  p)u  -  uoxx  -  vrxy  -  wtxz  -  k.Tx)  (3) 

G  =  (pv  ,  puv  -  ryx  ,  pv2  +  p  -  ayy  ,  pv w  -  ryz  , 

(e  +  p)v  -  uryx  -  voyy  -  wryz  -  K Ty)  (4) 

H  =  (pw  ,  pu W  -  Tzx  ,  pvw  -  Txy  ,  pw 2  +  P  -  <7ZZ  , 

(e  +  p)w  -  UTZX  -  VTzy  -  wazz  -  kTz )  .  (5) 

Metric  Coefficients 

The  metric  coefficients  of  the  generalized  coordinate  transformation  (1)  are  computed 
numerically  as  the  difference  interval  averages  of  the  mesh  paint  functions. 
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In  the  numerical  method  the  expressions  on  the  right  hand  side  of  equations  (6)  are  evalu¬ 
ated  by  central  differences  except  at  boundaries  where  one  sided  differences  to  the  interior 
are  employed.  The  Jacobian  of  the  transformation  J  is 

J  =x(c(yfJZ0  -  y^,zn)  +  yi{zr,x^  -  xnz 0) 

+  zt[x„yxi)  -x^)  (7a) 

=*«£*  +y«G  +  (7t) 

As  directly  interpreted  equations  (1),  (6)  and  (7a)  describe  a  finite  volume  method 
with  (6)  the  cell  face  area  and  (7a)  the  volume.  A  finite  (divided)  difference  upwind  method 
such  as  is  implemented  here  is  obtained  by  dividing  the  interval  averaged  metrics  by  one 
sided  difference  cell  volumes  as  given  by  equations  like  (7b)  for  each  coordinate  direction. 

Viscous  Fluxes 

The  viscous  terms  of  the  flux  vectors 

axx  =2 [n  +  fiT)(ux-  (lij  +  vu  +  wx)/3) 

Gyy  =  2(/l  +  t*T)(Vy  ~  (u  X  +Vy  +  W  x)  /  3) 

°zz  =  2(/i  +  Hr){ wx  -  (u*  +  Vy  +  wx)/3) 
rxu  =  Tvx  -  [n  +  Ar){u1(  +  v*) 

Tyx  =  Txy  -  (M  +  Hr){Vz  +  t Uy) 

Txz  =  Txx  =  (/X  +  Mr  )(«*  +  Ux) 

are  second-order  central  differences  in  space  and  first-order  implicit  in  time.  These  terms 
are  treated  with  the  Baldwin-Lomax11  thin  layer  approximation. 

When  turbulence  effects  are  present,  they  are  presently  modeled  with  the  well  known 
Baldwin  and  Lomax11  algebraic  mixing  length  model.  The  Baldwin-Lomax  algebraic 
model  defines  an  inner  and  outer  layer  formulation  for  the  turbulent  eddy  viscosity  in 
wall  bounded  shear  layers,  while  only  the  outer  layer  is  used  for  free  shear  layers  and 
wakes.  The  main  advantage  of  this  model  is  that  the  velocity  and  length  scales  are  easily 
defined  and  the  model  provides  reasonable  predictions.  This  model  has  been  modified 
following  the  ideas  of  Bradshaw12  for  near  wake  flows,  by  setting  a  lower  bound  value  for 
the  wake  law  eddy  viscosity  as  a  function  of  the  maximum  viscosity  value  across  the  wake. 
This  turbulence  model  has  been  used  effectively  to  simulate  2-D/axisymmetric  wakes  and 
jet  flows13,14,15.  The  flow  of  the  present  application  around  the  forebody  reentry  vehicle 
at  high  altitude  is  essentially  laminar.  The  turbulence  model  incorporated  in  the  code 
would  be  useful  in  the  simulation  of  the  wake  region  or  in  aerodynamics  problems  at  lower 
altitudes. 
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Inviscid  Fluxes 


The  inviscid  flux  terms  of  equation  (l)  are  approximated  with  CSCM  flux  difference 
eigenvector  splitting  in  an  extension  of  the  approach  of  reference  5.  First  order  upwind 
spatial  differencing  is  unconditionally  stable  with  backward  Euler  pseudo-time  iteration. 
For  higher  order  upwind  spatial  differences,  we  incorporate  flux  limiting  schemes  of  Yang, 
Lombard  and  Bardina19. 

The  CSCM-S  method6,7,  which  is  a  (globally)  single  data  level  relaxation  technique,  is 
employed  in  marching  along  the  ^-direction.  Within  a  £  constant  plane  at  each  marching 
step  the  implicit  equations  are  relaxed  using  the  pseudo  time  dependent  CSCM  method 
which  is  a  two  data  level  relaxation  technique  employing  DDADI  approximate  factorization 
in  the  xp  and  t)  directions.  The  resulting  implicit  method  of  planes  with  one  (additional) 
inner  iteration  at  each  spatial  (£)  marching  step  provides  a  faster  convergence  of  the 
solution  for  nonlinear  systems  of  equations.  While  the  processing  time  per  grid  point  per 
iteration  is  increased,  the  number  of  total  iterations  is  substantially  decreased;  especially 
when  the  ^-direction  is  aligned  with  the  main  flow  direction.  Of  equal  importance  as  we 
shall  discuss  in  a  later  section,  this  method  also  provides  a  substantial  improvement  in 
efficiency  of  data  management  of  large  3-D  data  bases  and  is  fully  vectorizable. 

Each  of  the  flux  difference  vectors  expressed  in  curvilinear  coordinates  is  split  ac¬ 
cording  to  the  direction  of  the  characteristic  wave  propagation  by  using  the  locally  one 
dimensional  similarity  transformation  that  diagonalizes  the  transformation  matrix  A  of 
the  flux  difference  vector  AF  by  the  conservative  variables  q,  through  the  relations 


AF  =  AAq  =  MTDT  1M~lAq 


(8) 


where  M  is  the  transformation  matrix  of  the  primitive  variable  differences  into  the  con¬ 
servative  variable  differences, 
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and  T  is  the  transformation  matrix  of  the  primitive  variable  differences  into  the  char- 


acteristic  variable  differences, 


where  the  upper  prime  denotes 

xn  ~  xn/ \j xn'2  +  Vn2  +  zo2> 

-Z/  -z.  /x2  5  ^2 

=  ijyts  +  tv  +  t, ; 

similarly,  in  all  the  other  variables.  The  volumetric  internal  energy  P  =  p/(7  —  l)  is  a 
primitive  variable  of  the  formulation. 

In  equation  (8),  to  effect  the  eigenvector  splitting  the  unit  matrix  is  split  into  three 
complementary  D+ ,  D~ ,  and  D°  diagonal  component  matrices.  The  diagonal  elements  of 
these  matrices  are  0’s  or  l’s  according  to  the  propagation  direction  of  the  characteristic 
waves,  as  given  by  the  sign  of  the  eigenvalues  of  the  transformation  matrix  A.  Thus 
D+  has  l’s  only  on  the  diagonal  coefficients  of  the  rows  corresponding  to  the  rows  of  A 
with  positive  eigenvalues;  D~  has  l’s  only  on  the  diagonal  coefficients  corresponding  to 
negative  eigenvalues,  while  D°  is  the  complementary  matrix  for  zero  eigenvalues.  This 
matrix  D°  which  is  not  separately  identified  in  other  numerical  methods  is  employed  to 
increase  robustness  and  stability.  The  eigenvalues  for  the  fluxes,  e.g.  in  the  ^-direction  arc 


W  ,  W  -f  ,  W  —  £c  ,  where  (9) 

w  =  (txpu  +  iypv  +  Zzpw)/p  (10) 

is  the  component  of  the  contravariant  averaged  velocity, 

€  =  (lx  +  +  h  )1/2  >and  (n) 

c  —  yp/'p  is  the  average  sound  speed.  (12) 


Consequently,  the  flux  difference  vector  of  equation  (8)  is  split  as 

A  F  =  A  F+  +  A  F~  =  (A  +  A  +  A~  A)q 


(13) 


where 

A±  =  ~MTD±T~lM- 1  (14a) 

and  A  is  the  forward- difference  operator.  Note  the  zero  eigenvalue  makes  no  contribution 
to  the  right  hand  side  flux  difference. 

By  analogy  to  the  well  known  result  for  the  scalar  advection  equation,  unconditionally 
stable  Euler  implicit  methods  are  constructed  of  the  flux  difference  split  pieces  by  upwind 
forward  differencing  with  pieces  associated  with  negative  eigenvalues  and,  similarly,  back¬ 
ward  differencing  pieces  with  positive  eigenvalues. 

To  improve  the  diagonal  dominance  we  impose  the  term 

A0  =MTD°T~1JI~1  (146) 

on  the  diagonal  of  the  left  hand  side.  The  method  is  well-posed  in  the  case  where  the  mag¬ 
nitude  of  any  eigenvalue  is  zero;  in  these  cases  the  approximation,  effectively  asserts  that 
the  characteristic  quantity  associated  with  a  zero  advective  eigenvalue  doesn’t  change17. 
The  effectively  frozen  eigenfunction  is  found  to  make  the  system  more  stable  and  produces 
smooth  solutions  particularly  through  strong  captured  shock  transitions. 

3-D  DDADI  Implicit  Diagonally  Dominant  Factorization 

For  computational  efficiency,  multi-dimensional  linearized  block  implicit  schemes  em¬ 
ploy  approximate  factorization  of  the  left-hand  side  matrices  to  permit  solving  coupled 
matrix  equations  in  one  family  of  coordinates  at  a  time. 

The  matrices  corresponding  to  zero  eigenvalues  and  that  are  added  to  I  on  the  left 
hand  side  are  shown  in  the  relationships  described  below.  Then,  the  3-D  unfactored  system 
sketched  as  a  first  order  method  on  the  right  hand  side  is 

{I  +  A°  +  L)6q  =  -Lq  (15) 

where  the  operators  L  and  A0  are  defined  as 

L  —  -f-  Lfj  +  (lGa) 

A°  =  A°e  +  A°  +  A%  (166) 

and 

L(  =  A%V€  +  AJ  A*  (17a) 

Ln  =  A+Vn  +  A~An  (176) 
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Ly,  —  A^Vy  +  Ay  (l"c) 

The  3-D  diagonal-dominant  approximate  factorization  DDADI5  of  the  unfactored  implicit 
upwind  matrix  equations  is  used  to  solve  the  system  of  equations.  The  resulting  relaxation 
scheme  for  the  method  of  planes  is  written  as: 

{I  +  A0  +  A(  + Ln  +  Ly)6q  = -Ln'n+1q  (18) 

where  the  absolute  value  matrix  operator  A £  on  the  diagonal  is  defined  as 

At  =  A+  -  A~  (19) 

The  missing  left  hand  side  contribution  of  the  off-diagonal  matrix  operator  —  A^  are 
operationally  included  in  the  right-hand-side  of  equation  (18)  in  the  updated  data  of  the 
symmetric  pair  of  space  marching  sweeps0.  Equation  (18)  is  solved  on  each  plane  in  the 
forward  and  backward  space  marching  sweeps  on  alternating  time  steps,  and  the  matrix 
operator  Z/n,n+1  indicates  that  the  data  belonging  to  all  previously  computed  planes  are 
already  updated. 

In  solving  equation  (18),  the  DDADI  approximate  factorization  of  the  left-hand-side  of 
equation  (18)  in  the  other  two  coordinate  directions  leads  to  the  following  block  tridiagonal 
equation  sequence: 

(D  +  Ln  -  An)6q*  =  -Ln'n+iq  (20) 

{D  +  Ly  -  Ay)Sq  =  -D8q*  (21) 

with 

qn+1  =  qn  +  6qn  (22) 

and  where  the  diagonal  matrix  operator  D  is  defined  as 

D  =  I  +  A0  +  A^  +  Aq  +  A^,  (23) 

Observe  that  the  absolute  values  of  the  eigenvalues  for  all  coordinate  directions  are  in  the 
diagonally  dominant  matrix  D,  equations  (20)  and  (21). 

The  Newton-Raphson  inner  iteration  is  applied  by  setting  up  the  equations,  solving 
and  updating  the  base  data  at  each  plane  twice  during  the  forward  and  backward  sweep 
procedure.  This  inner  iteration  is  central  to  the  improvement  in  fast  rate  of  convergence 
of  the  solution  of  the  nonlinear  system  of  equations.  The  nonlinearities  of  the  Navier- 
Stokes  equations  are  better  approximated  in  the  updated  locally  active  matrices,  and  the 
characteristic  waves  are  propagated  more  accurately6. 

Boundary  Conditions 

Implicit  characteristic  boundary  procedures  based  on  the  approach  of  reference  5 
are  applied  at  each  boundary  where  at  least  one  characteristic  wave  propagates  from  the 
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boundary  toward  the  interior  of  the  computational  domain.  The  boundary  approximations 
are  added  as  system  of  equations  similar  to  the  one  for  interior  grid  points  ol/equations 
(18);  except  for  the  absence  of  convection  matrices  A+  at  left  boundaries  and  A~  at  right 
boundaries  which  are  replaced  by  local  transformation  matrices,  similar  to  A0,  that  contain 
the  linear  representations  of  the  characteristic  variables  and  auxiliary  (imposed)  boundary 
conditions  in  terms  of  conservative  variables. 

For  supersonic  freestream  inflow  conditions  such  as  the  outer  boundary  of  reentry 
vehicles,  the  conservative  flow  variables  are  properly  frozen.  For  supersonic  outflow  condi¬ 
tions,  the  conservative  variables  are  properly  extrapolated,  i.e.  calculated  in  the  solution 
procedure,  from  the  interior  domain.  For  subsonic  outflow,  such  as  the  one  inside  of  the 
boundary  layer  zone,  invariance  in  the  static  pressure  is  invoked  by  the  solution  procedure 
differencing  toward  the  outflow  boundary;  however,  the  pressure  field  is  allowed  to  appro¬ 
priately  relax  since  the  solution  procedure  is  also  applied  differencing  along  the  outflow 
boundary.  For  viscous  wall  boundaries,  no-slip  conditions  are  applied  and/or  specified 
blowing/suction  mass  rates,  together  with  specified  wall  temperature  and/or  adiabatic 
wall  conditions.  At  symmetry  planes,  implicit  symmetric  approximations  of  the  flux  con¬ 
tributions  beyond  the  boundary  are  imposed  by  symmetry  relations  based  on  the  interior 
procedure  toward  the  boundary. 

Data  Storage  and  Management 

Three-dimensional  numerical  simulations  often  require  large  data  bases  of  the  order 
of  few  to  several  million  words  of  storage.  Due  to  the  single  data  level  relaxation  technique 
and  with  recomputed  metrics,  the  present  method  requires  only  eight  variables  per  grid 
point,  corresponding  to  the  three  coordinate  directions  and  the  five  dependent  variables. 
For  large  data  bases  on  the  CRAY-XMP  computer  that  we  have  been  using,  the  global 
arrays  are  stored  in  SSD  solid-state  devices.  These  devices  have  a  fast  transfer  rate  to 
the  central  memory  of  the  vector  processors  where  our  method  operates  on  only  a  few 
(five  for  the  second  order  method)  £  data  planes  at  a  time.  In  the  absence  of  the  single 
level  technique,  at  least  five  other  (residual)  variables  with  total  gnd  dimensions  will  be 
required.  Other  variables  may  also  be  stored  in  order  to  avoid  repetitive  and  unnecessary 
computation,  particularly  the  nine  metric  coefficients  and  one  Jacobian  of  the  coordinate 
transformation. 

With  first  order  differencing  in  the  ^-direction,  such  as  is  implemented  at  present,  the 
£  sweeps  require  loading  the  data  of  planes  J  —  1,  J ,  and  J  + 1  into  central  memory  in  order 
to  process  the  data  of  plane  J.  In  coding,  the  computation  of  the  flux  terms  and  matrix 
coefficients  are  fully  vectorized  in  the  £  plane.  The  boundary  conditions  and  solvers  are 
also  vectorized.  Once  the  data  of  plane  J  is  updated,  this  new  data  is  stored,  and  the  data 
of  plane  J  +  2  (on  a  forward  sweep)  is  loaded  into  centra!  memory  in  order  to  advance  the 
data  of  plane  J  +  1.  Thus,  the  complete  forward  £  sweep  requires  loading  the  data  base 
only  once  into  central  memory,  and  so  does  the  backward  £  sweep. 
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Appropriate  initialization  of  the  three-dimensional  data  is  required  in  order  to  obtain 
convergence  in  comparatively  few  numerical  iterations.  In  present  work,  solutions  were  first 
obtained  for  the  case  of  3-D  axisymmetric  flow  with  the  2-D/axisymmetric  flow  code5.  The 
3-D  solutions  are  obtained  within  a  few  hundred  iterations  with  about  9  10-5  cpu  seconds 
per  iteration  and  per  grid  point.  This  is  a  reasonable  fast  processing  time  considering  that 
this  is  DDADI  diagonal-dominant  approximate  factorization  scheme.  Then  the  solution  is 
run  with  the  3-D  code,  changing  only  the  free  stream  flowfield  conditions,  and  the  solution 
is  obtained  within  N  complete  alternating  sweeps  where  N  is  about  half  of  the  grid  points 
in  the  streamwise  direction.  The  three  dimensional  solutions  relaxed  in  about  60  global 
iterations  and  about  an  hour  of  CPU  time.  Simpler  initialization  procedures,  e.g.  based 
on  approximate  shock  fitting,  can  be  similarly  effective. 

Numerical  Results 

Numerical  results  of  the  hypersonic  flow  around  a  reentry  vehicle  (RV)  are  now  pre¬ 
sented.  The  geometry  of  the  vehicle  is  a  sphere-cone-cylinder-flare  body.  This  geometry 
has  previously  been  treated  as  Case  4  of  Kim  and  Lewis18.  Each  part  of  the  cone,  cylin¬ 
der,  and  flare  section  was  10  nose-radii  axial  lenght  with  a  half  angle  of  16°  in  the  cone 
surface  and  a  half  angle  of  12°  in  the  flare  surface.  The  free  stream  Mach  number  is  22, 
the  Reynolds  number  is  105,  Prandtl  number  is  0.7,  and  the  ratio  of  specific  heats  -y  =  1.4. 
A  cold  wall  boundary  condition  was  applied  with  Tw/Toa  —  7. 

Figures  la  and  lb  show  different  views  of  the  64x31  point  surface  mesh  on  the  reentry 
vehicle.  The  computational  mesh  has  64  points  in  the  streamwise  direction,  47  points  from 
the  body  wall  toward  the  outer  boundary,  and  31  points  in  the  circumferential  direction. 
The  mesh  has  been  constructed  by  using  efficient  2-D  algebraic  procedures22  based  on 
transfinite  interpolation  in  the  upper  and  lower  azimuthal  symmetry  planes,  shown  in 
Figure  lc,  and  a  similar  procedure  with  rotation  between  these  two  planes  in  order  to 
define  the  remaining  grid  planes,  as  is  shown  in  Figure  Id.  Figure  Id  shows  a  projection 
normal  to  the  axis  of  the  body  of  the  computational  mesh  plane  in  the  middle  of  the  flare 
section.  This  simple  procedure  is  very  effective  for  this  regular  geometry  to  define  the  grid 
with  minimum  distortions  and  gives  an  excellent  resolution  around  the  body  geometry 
and  corner  singularities.  The  grid  lines  are  normal  to  the  body  wall  and  follow  the  general 
geometry  of  the  bow  shock. 


Firstly,  Figures  2a,  2b,  and  2c  show  contour  plots  of  Mach  number,  normalized  pres¬ 
sure,  and  density  of  the  hypersonic  axisymmetric  flow  with  0°  of  angle  of  attack  along  a 
streamwise  plane.  The  pressure  is  normalized  with  its  free  stream  value.  The  bow  shock  is 
well  captured  and  the  contour  lines  in  the  shock  layer  are  smooth  and  free  of  oscillations 
and/or  instabilities.  The  expansion  waves  generated  at  the  cone-cylinder  corner  and  the 
weak  oblique  compression  shock  generated  in  the  cylinder-flare  corner  are  well  captured 
by  the  method  as  shown  in  the  pressure  contours.  The  non-alignment  of  the  grid  with  the 
bow  shock  is  reflected  as  some  small  changes  in  the  shock  distribution,  as  well  as,  in  the 
lower  resolution  in  the  flare  section.  The  use  of  adaptive  grid  techniques  should  produce 


optimal  results  in  this  kind  of  simulation. 

In  order  to  test  the  numerical  results,  we  show  in  Figures  3a,  3b,  and  3c  similar  contour 
plots  to  those  of  Figures  2,  but  obtained  with  the  2-D  code  for  3-D  axisymmetric  flows. 
The  agreement  is  excellent  between  the  flow  contours  of  the  3-D  and  the  axisymmetric 
3-D  simulations.  The  external  bow  shock,  the  expansion  fan,  and  recompression  shock  are 
almost  identical  in  both  simulations.  Minor  differences  can  be  seen  in  the  leeward  side  of 
the  cylinder-body  section. 

Figures  4a,  4b,  and  4c  show  the  contour  plots  of  Mach  number,  normalized  pressure, 
and  density  of  the  hypersonic  flow  at  10°  of  angle  of  attack  along  the  bisymmetry  plane. 
This  numerical  solution  was  obtained  after  running  60  global  iterations  or  complete  alter¬ 
nating  sweeps  over  the  axisymmetric  solution.  The  residuals  of  the  equation  of  motions 
decreased  at  least  3  orders  of  magnitude  and  the  solution  was  converged.  The  flow  struc¬ 
tures  are  well  captured  and  the  contour  lines  are  smooth  and  free  of  oscillations  and/or 
instabilities.  The  expansion  waves  generated  at  the  cone-cylinder  corner  and  the  weak 
oblique  compression  shock  generated  in  the  cylinder-flare  corner  are  well  captured  by  the 
method  as  shown  in  the  pressure  contours.  Also  evident  is  the  inflection  of  the  bow  shock 
with  the  interaction  of  the  expansion  fans  off  the  cone-cylinder  junction.  The  pressure  and 
density  contours  show  a  normal  section  of  the  recompresion  shock  off  the  cylinder-flare 
junction. 

Figure  5a  shows  the  projection  of  the  surface  body  mesh  onto  the  bisymmetry  plane, 
while  Figure  5b  show  the  contour  plot  of  the  normalized  pressure  on  this  plane.  The  effect 
of  the  free  stream  angle  of  attack  on  the  pressure  distribution  is  clearly  defined  in  this 
Figure  5b,  the  highest  contour  values  are  on  the  sphere  surface  in  the  windward  side  and 
the  smallest  contour  values  are  on  the  cylinder  section  in  the  leeward  side.  Figure  6a  shows 
the  surface  pressure  coefficient  versus  distance  along  the  body  surface.  The  horizontal  axis 
of  this  plot  is  centered  on  the  intersection  of  the  axis  of  symmetry  of  the  RV  with  the 
surface  of  the  sphere  section.  The  symbols  show  the  pressure  coefficients  at  10°  angle  of 
attack  along  the  body  surface  at  every  15°  from  the  bisymmetry  plane,  while  the  solid  line 
shows  the  pressure  coefficient  at  0°  angle  of  attack.  Figure  6b  shows  a  section  of  Figure  6a 
of  the  pressure  distribution  along  the  sphere-body  section.  The  effects  of  the  free  stream 
angle  of  attack  is  well  captured.  There  is  some  lack  of  smoothness  near  the  pole  of  the  grid 
topology  due  to  an  inconsistency  of  the  grid  spacing  with  the  cross  pole  difference  scheme 
that  does  not  solve  at  the  pole.  The  pole  should  have  been  “buried”  at  half  cell  spacing 
in  the  mesh. 

Figures  7  show  normalized  pressure  contours  of  cross  sectional  mesh  planes  projected 
onto  the  normal  plane  to  the  axis  of  body  symmetry.  Figures  7a,  7b,  and  7c  correspond  to 
midplane  section  of  the  cone,  cylinder,  and  flare  sections,  respectively,  at  10°  of  angle  of 
attack.  On  the  other  hand,  Figure  7d  correspond  to  midplane  section  of  the  flare  section 
at  0°  of  angle  of  attack.  Similarly,  Figures  8  show  the  density  contours  corrsponding  to  the 
respective  Figures  7.  These  figures  show  the  three-dimensionality  of  the  flow  structures, 
the  development  of  the  bow  shock,  the  3-D  zone  of  expansion  fan,  and  3-D  zone  of  recom- 
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pression  shock.  The  complexity  of  the  flow  is  shown  to  be  well  captured  by  the  CSCM3D 
numerical  method. 

Concluding  Remarks 

The  3-D  CSCM  upwind  compressible  Navier-Stokes  algorithm  has  been  improved  and 
extended  to  simulate  external  hypersonic  flows  at  angle  of  attack.  The  method  combines 
the  best  features  of  storage  and  processing  efficiency  of  PNS  space  marching  procedures 
with  the  generality  of  time  dependent  techniques  to  solve  flows  with  elliptic  and  streamvvise 
separated  zones.  The  3-D  DDADI  diagonal-dominant  approximate  factorization,  together 
with  well  posed  implicit  boundary  approximations,  and  the  frozen  eigenfunction  approxi¬ 
mation  at  zero  eigenvalue  crossings  provide  robust  stability  and  fast  relaxation  to  steady 
state.  Numerical  results  for  Mach  22  flow  at  0°  and  10°  angle  of  attack  of  a  sphere-cone- 
cylinder-flare  reentry  vehicle  are  presented.  The  numerical  results  are  consistent  with 
axisymmetric  simulations  and  with  previous  experience  in  external  hypersonic  flows. 
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Figure  lc.  Computational  mesh  on  bisymmetry 
plane  of  RV  veh icl e . 


Figure  2c.  Density  contour  lines 
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Figure  5a.  Surface  mesh  body  projected  onto 
bi symmetry  plane. 


Figure  5b.  Normalized  pressure  contour  lines 
of  hypersonic  flow  at  10°  angle  of  attack, 
correspond i ng  to  figure  5a. 


-no-Joo  -no  -no  -no  -too  -so  oo  so  wo  no  joo  »o  soo  jso 
VKn 

Figure  6a.  Surface  pressure  coefficient  along  RV 
body.  Symbols  show  predictions  along  planes  at 
0M5‘\30\45<\6O°,75'\and  90°  off  the  bisym- 
metry  plane  of  hypersonic  flow  at  10°  angle  of 
attack.  Solid  line  shows  prediction  of  axi- 
symmetric  flow  at  0°  angle  of  attack. 


Figure  6b.  Surface  pressure  coefficient  along 
the  sphere-body  section.  Symbols  correspond 
to  figure  6a. 
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8.  Interactions 

The  research  described  in  this  report  has  been  partially  presented  in  seven  technical 
meetings  in  the  U.S.  and  abroad.  The  research  in  simplified  algebraic  grid  generation  is 
partially  described  in  Section  1  reference  2  of  this  annual  report. 

The  solution  adaptive  grid  work  on  multiple  patch  meshes  has  been  given  in  a  paper, 
Accurate  Numerical  Simulation  of  Supersonic  Jet  Exhaust  Flow  with  CSCM  on  Adaptive 
Overlapping  Grids,  AIAA-87-0465,  AIAA  25th  Aerospace  Sciences  Meeting,  Reno,  NV, 
January  12-15,  1987.  The  work  has  also  been  presented  in  a  talk  by  C.K.  Lombard  at  the 
Second  Nobeyama  Workshop  on  Supercomputering  in  Fluid  Dynamics,  Nobeyama,  Japan, 
September  7-10,  1987. 

Work  on  the  CSCM  upwind  method  expressed  as  a  finite  volume  method  on  node 
bounded  cells  in  staggard  grids  has  been  partially  presented  in  a  paper,  Accurate,  Efficient 
and  Productive  Methodology  for  Solving  Turbulent  Viscous  Flows  in  Complex  Geometry, 
for  the  10th  International  Conference  on  Numerical  Methods  in  Fluid  Dynamics,  Beijing, 
China,  June  23-27,  1986. 

Work  on  nonlinearly  stable  high  order  biased  upwind  methods  has  been  partially  pre¬ 
sented  in  a  paper,  Uniformly  Second  Order  Accurate  ENO  Schemes  for  the  Euler  Equations 
of  Gas  Dynamics,  AIAA-87-1166-CP,  AIAA  8th  Computational  Fluid  Dynamics  Confer¬ 
ence,  Honolulu,  Hawaii,  June  9-11,  1987. 

The  CSCM  shock  capturing  scheme  with  improved  diagonal  dominance  through  frozen 
eigenfunctions  for  zero  eigenvalue  crossings  has  been  presented  with  application  of  3-D 
hypersonic  blunt  body  flow  in  the  paper,  Three  Dimensional  Hypersonic  Flow  Simulations 
with  the  CSCM  Implicit  Upwind  Navier-Stokes  Method,  AIAA-87-1114-CP,  AIAA  8th 
Computational  Fluid  Dynamics  Conference,  Honolulu,  Hawaii,  June  9-11,  1987. 

Recent  application  of  the  combined  algebraic  grid  generation  and  surface  geometry 
definition  tools  with  the  segmented  patched  2-D/axisymmetric  CSCM  Navier-Stokes  solver 
(developed  under  AFOSR  support)  to  the  design  modification  of  the  throat  region  of  the 
Mach  14  nozzle  of  the  NASA-Ames  3.5  ft  hypersonic  tunnel  are  reported  in  the  paper, 
Aerodynamic  Design  Modification  of  a  Hypersonic  Wind  Tunnel  Nozzle  by  CSCM  with 
High  Order  Accuracy,  AIAA-87-1896,  AIAA/SAE/ASME/ASEE  23rd  Joint  Propulsion 
Conference,  San  Diego,  CA,  June  29-July  2,  1987. 

9.  New  Discoveries 

The  CSCM  method  on  node  bounded  finite  volume  computational  cells  in  a  staggard 
grid  is  a  new  class  of  flux  difference  split  upwind  schemes. 


